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Abstract 

The  United  States  Air  Force,  Department  of  Defense  and  commercial  industry 
have  recognized  the  great  value  of  near-earth  space  development,  specifically  in  satellites 
for  use  in  communications,  ground  and  space  surveillance  and  more  active  roles. 
However,  resolution,  or  the  primary  optic’s  diameter,  has  been  a  limitation,  especially  for 
ground  surveillance.  Deployable  optics  has  been  investigated  to  allow  larger  optics  in 
space  and  membrane  optics  has  received  increasing  attention  recently.  The  membrane’s 
flexible  nature  requires  some  passive  and  possibly  active  control  to  reduce  optical 
distortion  caused  by  manufacturing,  deployment,  or  other  effects  during  use. 

Piezoelectric  surface  controllers  are  one  option  to  actively  control  the  membrane  on  the 
order  of  optical  measurements  (micron  displacement  or  less).  Multiple  configurations  of 
transverse  displacements  are  feasible  depending  on  the  piezo  zone  locations  and 
activation. 

The  current  thrust  of  industry  is  reducing  the  effort,  time  and  cost  of 
manufacturing  and  testing  through  use  of  computerized  modeling  and  simulation; 
therefore,  this  was  investigated  for  a  membrane  mirror  and  piezoelectric  combination. 
Prior  experiments  using  6-  inch  diameter  membranes  have  been  conducted  with  an 
axisymmetric  piezoelectric  material  layer  on  the  non-optical  surface.  Various  voltage 
differentials  were  applied  to  the  piezo  and  the  transverse  displacement  was  measured.  A 
finite  element  code,  using  perturbation  techniques,  was  written  in  MATLAB  and  tested  to 
check  the  feasibility  of  using  computer  models  for  the  micro- displacements  occurring 
with  the  membrane- piezo  lay-up.  The  computer  program  considered  was  developed  for 


IV 


axisymmetric  conditions;  however,  in  many  cases,  these  conditions  tended  to  dominate. 
Under  these  conditions,  the  finite  element  code  produces  results  that  represent  the 
axisymmetrically  reduced  experimental  data. 
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AXISYMMETRIC  OPTICAL  MEMBRANE  MODELING  BASED  ON 
EXPERIMENTAL  RESULTS 

I:  Introduction  and  Background 

Space  Optics  And  Other  Uses 

The  United  States  Air  Force  and  Department  of  Defense,  as  well  as  other 
government  and  civilian  organizations,  have  recognized  the  need  for  larger  space 
telescopes  for  investigating  the  far  reaches  of  the  universe  and  discovering  parts  of  the 
Earth  otherwise  impossible  to  see  (Marker/Jenkins,  2001).  The  motivation  behind  this  is 
to  avoid  atmospheric  disturbances  when  searching  the  stars  and  for  the  coverage  the 
altitude  gives  when  viewing  the  Earth.  However,  the  resolution  is  reduced  due  to  this 
altitude.  As  a  general  rule,  the  resolution  is  directly  related  to  the  diameter  of  the  primary 
mirror.  (Sobers,  2002)  Solid  optics  has  been  the  chief  type  used  on  these  systems  in  the 
past.  Due  to  limitations  in  payload  space  and  weight,  these  optics  are  not  feasible  over  a 
few  meters  in  diameter. 

The  military  community  and  DoD  are  most  interested  in  an  earth  observing  ability. 
Past  technologies  to  accomplish  this  have  focused  on  solid  optics.  The  most  familiar  and 
well-known  optical  quality  telescope  is  the  Hubble  (possessing  a  2.4-  m  solid  primary), 
used  to  photograph  the  stars  and  distant  galaxies.  If  this  were  turned  toward  the  ground  it 
would  have  resolution  capabilities  of  0.15  m.  (AFRL/DEBS  TMR,  2002).  There  are  also 
other  satellite  platforms,  one  such  named  IKONOS,  capable  of  ground  resolution  of  0.6 
m.  A  25-  m  membrane  placed  at  12000  km  could  potentially  have  a  ground  resolution  of 
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0.3  m,  while  having  a  swath  width  of  840  km  verses  the  Hubble  at  only  about  1.2  m. 

These  calculations  are  completed  using  the  following  equation: 

Resolution=  1.22XR/D  (1) 

where  X  =0.5  pm,  D  is  the  mirror  diameter  and  R  is  the  range  of  the  system. 

While  serving  some  of  the  needs  for  which  they  were  designed,  a  higher 
resolution  system  is  desired.  Achieving  this  requires  redesigning  the  primary  mirror  to 
better  suit  the  needs  and  capabilities  of  modern  space  systems.  Some  ideas  for  the 
redesign  include  a  segmented  solid  or  a  membrane  optic.  Both  have  their  own  limitations 
and  strengths,  and  both  are  being  developed  at  some  level  in  and  out  of  the  DoD; 
however,  only  one  is  investigated  in  this  study.  (Magee,  Carlin) 

Segmented  solid  optics  have  been  investigated  as  a  solution  for  quite  some  time 
as  a  way  to  increase  the  full  primary  mirror  size,  but  have  the  same  or  less  payload  size. 
This  requires  deployment  of  the  optic  in  some  fashion,  such  as  a  flower  opening  or 
fanning  and  positioning  of  the  segments  to  within  fractions  of  a  wavelength.  They  offer 
the  opportunity  to  increase  the  optic  size  using  current  solid  optics  technology  mated  with 
current  control  systems  to  bring  the  segments  into  the  proper  positioning.  (Li,  1998) 
Segmented  optics  share  some  of  the  same  shortcomings  as  solid  optics,  including,  but  not 
limited  to  the  weight  of  all  the  segmented  pieces,  the  complications  of  deployment  when 
on  station  and  the  expensive  and  timely  procedure  of  polishing  the  mirror  pieces. 

(Magee,  2002) 

In  recent  years,  there  have  been  space  telescope  experiments  using  polymers;  one 
such  was  NASA’s  Inflatable  Antenna  Experiment  (IAE),  which  is  a  pressurized  lenticular 
about  14  meters  in  diameter.  It  was  designed  as  a  radio  antenna  and  its  main  purpose  was 
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verification  of  the  deployment  mechanics.  (L’Garde,  2004)  Figure  1  shows  the  inflated 
antenna  in  the  laboratory.  This  technique  has  little  application  in  the  visible  and  shorter 
wavelengths  since  the  front  canopy  causes  diffraction  of  the  light  wave  before  reflecting 
off  the  primary,  which  itself  is  not  of  diffraction  limited  quality;  however,  the  success  of 
IAE  has  lead  to  the  exploration  of  polymers  as  a  probable  primary  mirror. 


Figure  1.  Laboratory  IAE 


Following  on  project  successes  such  as  the  IAE  and  the  continued  desire  to 
advance  current  capabilities,  thin  polymers  have  been  researched  in  greater  depth  to 
discover  their  possible  uses  in  optical  systems.  There  have  been  many  speculations  on 
the  end  product  of  the  membrane  mirror.  For  instance,  reconnaissance  satellites,  relay 
satellites,  and  ground-based  astronomy.  (Rotge,  2002)  Recently,  there  has  even  been 
discussion  on  the  use  of  the  membranes  in  a  Relay  Mirror  Airship  used  in  conjunction 
with  the  Airborne  Laser  (ABL)  to  increase  the  coverage  and  decrease  the  absolute 
number  of  ABL  platforms  that  need  to  be  airborne  at  any  one  time.  Lor  these  reasons  and 
to  help  make  the  previous  possible,  membrane  optics  are  selected  to  be  further  studied  for 
large  aperture  space  telescope  development. 
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Membrane  Optical  Considerations 


Membrane  optics  provide  many  advantages  over  traditional  solid  optics.  First, 
they  are  usually  cheaper  and  more  time  efficient  to  manufacture,  since  there  is  only  one 
manufacturing  stage,  eliminating  the  polishing  process.  (Rotge,  1999)  They  have  the 
opportunity  to  reduce  the  payload  weight  by  many  magnitudes,  since  their  thickness  is 
essentially  nil  compared  to  their  diameter,  reducing  the  cost  of  the  launch.  Lastly,  the 
membrane  has  the  opportunity  to  be  adapted  to  many  other  missions  with  greater  ease 
than  a  solid  optic.  One  example  is  the  use  of  the  membrane  as  an  extremely  thin  optical 
window,  thus  preventing  the  light  from  slowing  or  being  diffracted  in  thicker  glass 
windows.  Of  course,  membrane  optics  share  deployment  complications  with  their 
segmented  cousins,  as  well  as  their  own  limitations,  such  as  surface  deflections  caused  by 
slewing  the  structure  and  manufacturing  the  membrane  to  near  optical  tolerances  over 
multi- meter  sections  without  polishing;  however,  the  benefits  outweigh  the  difficulties 
and  many  of  the  limitations  have  been  addressed  and  some  of  them  are  already  overcome. 
(Marker,  2002)  The  chief  reason  for  investigating  membrane  optics  is  overall  weight 
reduction  of  the  total  system  by  using  a  membrane  whose  weight  in  comparison  to  the 
support  structure  is  nil;  thereby,  reducing  launch  cost  and  increasing  the  opportunity  to 
place  into  a  higher  orbit. 

Before  membranes  could  be  considered  further,  better  manufacturing  techniques 
were  needed  to  produce  the  optical  tolerances  without  polishing.  A  few  techniques  were 
investigated  and  one  method  commonly  used  is  a  spin  casting  method  that  is  currently 
able  to  manufacture  flat  membrane  mirrors  on  the  order  of  1-m  diameter  with  X  /20 
diffraction  over  the  whole  span,  meaning  that  there  is  a  thickness  difference  of  about  30 
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nm  across  the  membrane  for  a  green  light  laser  (A,  =  633  nm).(SRS  website)  The 
material  science  development  of  these  membranes  (CP-1)  was  completed  through  USAF 
contracts  at  the  AFRL  Directed  Energy  Directorate  and  NASA.  This  work  is  ongoing  to 
accomplish  this  same  feat  on  the  multi- meter  level.  (Rotge,  2000) 

The  manufacturing  of  any  optic  is  usually  developed  around  the  global  shape 
while  polishing  the  glass  or  metal  will  reduce  the  amount  of  local  hills  and  valleys.  This 
affects  the  path  of  the  incident  light  on  the  optic.  As  collimated  light  encounters  a 
surface,  it  behaves  in  two  fashions.  First,  there  is  a  global  reflection  and  refraction  that 
behaves  geometrically,  i.e.  if  the  light  is  incident  to  the  surface  at  a  certain  angle,  it  will 
reflect  (bounce  off  surface)  at  the  same  angle  and  refract  (enter  the  medium)  at  an  angle 
based  on  Snell’s  Law: 

nj  sin  0(.  =  nt  sinO,  (2) 

where  nLt  is  the  index  of  refraction  (Hecht) 

The  reflection  is  of  more  concern,  since  any  light  that  enters  the  membrane  will 
mostly  exit  the  other  side  and  a  thin  layer  of  a  specific  wavelength  reflective  coating  can 
prevent  much  of  the  refraction.  The  second  reflection  characteristic  requires  looking  at 
the  more  localized  surface,  since  in  optics  the  light  is  usually  a  specific  laser  wavelength, 
usually  a  shorter  one,  such  as  green  (633  nm)  in  this  case.  The  best  way  to  explain  what 
happens  locally  is  to  examine  Figure  2. 
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Figure  2.  Local  Reflection 

As  it  can  be  seen,  the  hills  or  valleys  will  cause  a  portion  of  the  light  to  reflect  at  a 
different  angle  from  its  incident  angle.  This  can  cause  quite  a  bit  of  loss  of  photons, 
thereby  reducing  the  light  entering  into  the  rest  of  an  optical  system.  However,  there  is 
a  way  to  prevent  or  correct  for  this  phenomenon.  One  current  method  is  real-time 
holography,  a  method  that  has  been  proven  for  many  wavelengths  of  diffraction  off  a 
membrane  (Gruneisen,  1999).  The  basic  explanation  of  this  method  is  to  use  two  beams 
of  light,  one  simulating  a  point  source  at  infinity  and  another  beam  that  is  “separated” 
into  a  few  sets  of  beams.  The  point  source  and  part  of  the  second  beam  are  collimated 
and  reflected  off  the  aberrated  membrane.  The  point  source  beam,  post  reflection,  is 
interfered  with  a  piece  of  the  original  second  beam  and  combined  with  another  piece  of 
the  original  second  beam  at  a  Spatial  Light  Modulator  (SLM),  where  the  reflected  second 
beam  meets  them.  The  SLM  is  where  the  correction  is  completed  and  the  light  is 
reflected  back  to  a  camera.  This  setup  can  be  seen  in  Figure  3  and  more  information  is 
supplied  in  the  literature  (Gruneisen,  1999). 
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aberrated  membrane 


Figure  3.  Holographic  Setup 


The  current  holography  methods  have  the  capability  to  correct  for  many 
wavelengths  of  diffraction  from  any  mechanism,  from  physical  distortion  of  the 
membrane  due  to  deployment  to  manufacturing  deficiencies.  There  is  also  the  ability  to 
correct  for  dynamic  changes  coming  about  from  slewing  or  other  dynamic  changes. 
Although  the  holography  can  correct  many  wavelengths,  it  was  recognized  early  that  this 
may  still  not  be  enough  and  an  actual  mechanical  correction  may  be  needed  as  well. 
While  there  have  been  many  ideas,  including  full  actuator  sets  as  was  planned  for  use 
with  segmented  optics  and  MicroElectroMechanical  Systems  (MEMS),  another  relatively 
inexpensive  and  seemingly  non- intrusive  idea  is  the  use  of  Polyvinylidene  fluoride 
(PVDF)  materials,  which  may  be  capable  of  active  control  through  a  voltage  differential 
placed  across  a  very  specific  surface  area  on  the  membrane.  Multiple  zones  can  be 
etched  out  and  therefore  be  activated  separately  or  in  conjunction  to  reshape  the 
membrane  to  suit  the  needs  of  the  mission,  the  manufacturing  deficiencies  or  for  any 
other  reason  that  the  membrane  needs  to  be  tweaked. 
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Currently,  most  membrane  manufacturing  techniques  have  concentrated  on 
developing  optical  quality  flats,  thin  films  with  no  curvature  from  the  manufacturing 
process.  (Rotge,  1999)  Within  the  earth  environment,  a  vacuum  could  be  applied  behind 
the  membrane  to  pull  it  into  a  convex  shape,  but  this  idea  has  two  problems,  it  takes  on  a 
spherical  shape  rather  than  a  parabolic  shape  usually  desired  and  it  is  not  feasible  in  the 
space  vacuum.  The  solution  is  a  mirror  that  is  globally  manufactured  to  the  desired  net- 
shape,  a  mirror  that  already  has  the  global  convex,  parabolic  shape  immediately  after 
manufacturing.  This  mirror  would  also  require  the  same  optical  tolerances  at  the  small- 
scale  to  avoid  diffraction  of  the  incoming  light.  Fortunately,  the  technology  used  for  flat 
membrane  manufacturing  is  applicable  to  net- shape  membranes.  (Marker,  2002)  The 
main  difference  is  the  casting  mold  shape,  namely  flat  or  a  specific  net- shape.  An 
advantage  to  this  technique  is  the  elimination  of  an  optically  toleranced  mold,  since  the 
air  surface  is  the  optical  surface. 

These  net- shape  membranes  provide  a  versatile  range  of  uses  from  space-based 
primary  mirrors  to  various  ground  uses.  One  concept  currently  in  development  is  to  use 
the  membrane  as  an  optical  coating  on  ceramic  mirrors,  thereby  reducing  the  need  to 
polish  the  mirror.  This  in  effect  reduces  time  from  start  of  mirror  manufacturing  to  “first 
light”,  or  when  the  mirror  can  first  be  illuminated  in  use.  As  with  many  other  types  of 
solid  mirrors,  the  difficulty  and  costly  procedure  of  polishing  the  ceramic  mirror,  a 
materially  hard  substance,  to  optical  tolerances  is  eliminated.  These  dual  layer  mirrors 
can  be  used  in  ground  applications  and  would  provide  the  capability  to  use  the  bottom 
layer  (metal  or  ceramic)  as  the  mold. 
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Another  challenge  of  the  membrane  is  the  lack  of  rigidity;  therefore,  a  Middle  or 
Geosynchronous  Earth  Orbit  (MEO/GEO)  is  optimal  to  avoid  atmospheric  disturbances, 
as  well  as  other  concerns  at  a  Low  Earth  Orbit  (LEO).  From  these  altitudes,  a  much 
larger  area  of  the  world  could  be  covered  by  a  single  satellite,  although  larger  primaries 
would  have  to  be  manufactured  to  meet  resolution  requirements.  If  local  defects  are  kept 
to  A.  /20  optical  tolerances  over  3-m  areas  of  the  membrane,  a  30- m  membrane  mirror 
could  potentially  be  feasible  for  the  visible  wavelengths;  especially  with  optical 
correction  methods  that  are  available  and  being  developed. 

Membrane  optics  have  many  unique  qualities  over  traditional  glass,  metal  or  other 
solid  optics.  The  first  one  that  normally  comes  to  mind  is  the  small  thickness-to-diameter 
ratio.  In  most  traditional  optics,  the  thickness-to-diameter  ratio  has  been  figured  at  about 
1:6  to  ensure  structural  and  dimensional  stability.  (Carlin)  Weight  reduction  of  solid 
optics  have  been  completed  through  use  of  honeycombing  or  bonding  a  face  sheet  to  a 
back  plate;  however,  there  are  inherent  difficulties  with  these  methods,  such  as  the 
fragility  of  the  mirror  and  out  of  plane  bending  after  honeycombing.  For  this  reason 
another  method  is  desired. 

A  membrane  optic,  on  the  other  hand,  can  be  manufactured  through  a  curing 
process  that  can  be  completed  on  the  order  of  hours  or  days  rather  than  weeks  or  months, 
to  the  required  optical  tolerances  at  thickness-to-diameter  ratios  of  1:100000  for  a  1-m 
membrane.  This  involves  an  extremely  stable  manufacturing  system  to  achieve  a  X  /20 
Root- Mean-Square  (RMS)  geometrical  figure,  the  surface  roughness  difference  across  a 
specific  portion  of  the  total  membrane.  This  may  be  best  described  using  a  sine  wave.  A 
perfect  sin(  0  )  wave  has  a  Peak-to-Valley  value  of  2  X  and  an  RMS  value  of  .707  X  . 
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Imagine  a  sine  wave  with  some  positive  error  at  the  first  positive  peak.  This  would  drive 
the  P-V  to  increase  the  total  amount  of  the  error  and  the  RMS  value  to  increase  by  an 
average  value  across  as  defined  by  the  definition  of  the  calculation.  Notice,  the  RMS 
value  describes  the  average  of  the  local  wavefront  in  comparison  to  the  global  wavefront, 
rather  than  describing  the  maximum  local  defects  like  the  Peak- to- Valley,  which  is  the 
difference  between  the  tallest  peak  and  the  deepest  valley. 

It  would  be  ideal  to  have  these  specifications  across  the  whole  membrane  at  large 
diameters,  but  is  unlikely  to  occur,  whether  for  membranes  or  glass.  To  envision  this, 
imagine  a  mirror  designed  for  green  light  reflection  (633  nm);  the  X  /20  requirement 
would  mean  an  average  surface  difference  nearly  30  nm  across  the  whole  front  surface. 

To  further  emphasize  this,  this  corresponds  to  a  ratio  of  surface  roughness  to  10- m 
diameter  mirror  of  30  x  10" 10  m  or  approximately  l/8000th  the  diameter  of  the  human  hair! 
Obviously,  the  larger  the  mirror  the  more  difficult  this  task  is  to  complete;  however,  if 
this  requirement  can  be  accomplished  over  any  3-m  diameter  area  of  the  membrane,  other 
correction  methods  (optical  holography  and/or  active  control)  may  supply  enough 
correction  to  meet  many  optical  missions. 

The  focus  of  many  projects  has  been  modeling  and  simulation  to  save  cost  and 
time  from  manufacturing  and  to  test  new  ideas  based  on  a  prior  idea.  This  is  in  the  early 
stages  for  membrane  mirrors  to  discover  how  they  will  respond  to  specific  forces  or 
control  mechanisms.  There  is  a  crucial  difference  between  finite  element  modeling  of 
membranes  and  that  of  other  structures,  namely,  the  order  of  deflections.  Deflections  in 
most  structures  are  usually  on  order  of  at  least  100s  of  microns  (0.1  mm),  if  not  much 
more;  however,  due  to  the  precision  of  an  optical  system,  a  finite  element  code  that  can 
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handle  nanometer  level  displacements  is  required.  In  earlier  completed  experimental 
work,  the  deflection  of  a  *4”  thick  membrane  instigated  by  a  PVDF  layer  was  on  the  order 
of  2  X  or  about  a  micron.  (Sobers,  2002)  This  is  usually  the  high  end  displacements  and 
fractions  of  one  wave  deflection  may  be  necessary  within  the  code.  Current  commercial 
codes  don’t  possess  this  capability  or  all  the  mechanics  of  the  PVDF  layer;  therefore,  the 
theory  was  researched  and  a  finite  element  code  was  developed  based  on  this  theory. 
(Rogers,  2001) 

The  scope  of  this  thesis  is  a  culmination  of  one  graduate  thesis  and  one 
dissertation.  The  thesis  dealt  with  experimental  results  of  activating  a  PVDF  layer  and 
how  it  affected  the  shape  of  the  membrane.  The  dissertation  involved  the  development  of 
a  perturbation  technique  finite  element  code.  This  study  covers  the  use  of  this  code  to 
check  it  against  the  reality  of  the  experimental  results.  This  work  will  be  covered  in  a 
general  sense  in  the  next  two  chapters. 
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II:  EXPERIMENTAL  DATA 


Overview 

As  stated  before,  all  experimental  testing  was  completed  under  a  separate  thesis. 
The  goal  was  to  improve  the  control  capability  of  an  optical  surface.  Past  systems  were 
limited  to  global  shape  control  or  small,  expensive  microelectricalmechanical  (MEM) 
devices.  To  improve  on  these  techniques,  four  mirrors  were  manufactured  using  a 
polymer-based  optical  surface  in  conjunction  with  two  types  of  bonded  materials.  The 
first  was  constructed  of  a  piezo-ceramic  bonded  to  a  copper-clad  circuit  board.  The  one 
investigated  in  this  study  was  a  stretched  PVDF  membrane  bonded  to  an  aluminum  ring. 
The  polymer  is  a  two  part  liquid  silicone  that  was  cured  in  the  ring  on  top  of  the  PVDF 
(see  Figure  4). 


Figure  4.  Cured  Polymer  in  Mount 

The  mirrors  were  tested  using  a  Shack- Hartmann  wavefront  sensor  to  determine  surface 
flatness.  The  surface  flatness  measurements  were  used  to  test  the  system  and  verify  it 


13 


was  setup  properly.  They  are  used  to  normalize  the  deflection  plots,  since  capturing  the 
deflection  and  not  the  initial  manufacturing  wavefront  is  desired. 

Mirror  Construction 

The  membranes  were  constructed  of  two  materials,  namely  GE  Silicone 
RTV615  and  a  piece  of  PVDF.  They  were  given  the  designation  Ml  and  M2.  Both 
mirrors  used  the  same  manufacturing  technique  with  a  slight  difference  in  the  PVDF. 

The  stretching  system  consisted  of  an  aluminum  ring  and  rubber  O-ring,  aluminum 
faceplate  and  four  bar- clamps. 

The  membrane  manufacturing  started  with  the  PVDF  being  placed  between  the  o- 
ring  and  faceplate.  The  bar-clamps  were  placed  over  the  outside  edge  of  the  faceplate 
and  tightened  until  the  PVDF  was  taut.  Epoxy  was  applied  to  the  aluminum  mounting¬ 
ring  on  the  middle  surface  between  the  two  grooves  and  then  place  on  the  stretched 
PVDF.  To  further  ensure  an  even  and  stretched  membrane,  weights  were  place  on  top  of 
the  mounting  ring.  The  PVDF  was  later  cut  out  from  the  stretching  system  after  the 
epoxy  cured. 

The  initial  purpose  of  Ml  was  to  test  the  mounting  techniques;  however,  the 
etching  procedures  were  later  developed  and  a  control  pattern  was  etched  on  the  back 
surface  after  mounting.  This  process  was  reversed  in  M2.  The  etched  pattern  was 
completed  before  the  PVDF  was  mounted.  Figure  5  shows  the  actual  M2  PVDF  layer 
and  Figure  6  shows  a  layout  of  the  activation  zones.  Notice  the  leads  to  each  of  the  zones. 
Zone  7  is  of  most  interest  due  to  its  perfectly  axisymmetric  location. 
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Figure  5.  Etched  PVDF  film 


Figure  6.  M2  Activation  Zones 


To  learn  more  about  the  etching  procedure,  refer  to  Smart  Structures  for  Control  of 
Optical  Surfaces  Appendix  A  (Sobers,  2002). 

In  both  cases  after  the  mounting  procedure  was  complete,  the  mirror  was  spray- 
painted,  silicone  rubber  primer  applied  and  a  layer  of  GE  Silicones  RTV615  was  poured 
into  the  mounting  ring  on  the  opposite  side  of  the  electrodes.  A  second  layer  of  RTV  was 
poured  due  to  curvature  that  developed  during  curing.  The  combination  of  the  two 
polymer  layers  was  measured  to  be  approximately  6  mm  Ml  and  4  mm  for  M2.  The 
PVDF  for  M2  was  slightly  thicker  at  5 2- pm  verses  32-pm  for  Ml. 
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Test  Methodology 

The  surface  flatness  measurements  were  completed  using  a  20  mW  helium- neon 
laser  with  a  wavelength  of  633  nm  and  a  WaveScope®  Shack- Hartmann  wavefront 
sensor  built  by  Adaptive  Optics  Associates.  The  data  outputted  was  used  to  calculate  a 
surface  plot,  synthetic  interferometric  fringe  pattern  and  the  first  35  coefficients  of  the 
Zernike  polynomial  set.  The  surface  data  and  Zernike  polynomials  were  exported  to  be 
plotted  in  MATLAB®  to  enhance  the  display. 

Using  various  optical  path  setups  and  neutral  density  filter  ranges,  the  surface 
flatness  was  measured  for  each  mirror.  The  detail  of  this  method  can  be  found  in  Chapter 
3  (Sobers,  2002).  After  the  surface  flatness  was  measured  and  beam  expansion  rates 
were  determined,  the  actual  testing  of  the  membrane  mirror  was  completed.  The  Shack- 
Hartmann  was  used  to  measure  the  test  and  reference  surfaces  separately.  This  was 
accomplished  in  this  manner  since  the  reference  surface  has  a  much  greater  reflection 
than  the  test  specimen.  It  also  removes  the  majority  of  test  equipment  bias  and  presents 
an  opportunity  to  compare  different  specimens.  Various  figures  of  the  setups  could  be 
seen  in  Figure  7 -Figure  8. 
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Figure  7.  Experimental  Filter  Setup 


Figure  8.  Beam  Expansion 

Along  with  different  setups,  a  few  different  types  of  equipment  were  used.  The 
Monolithic  Lenslet  Modules  MLM  lenslets  size  used  determines  the  fidelity  of  the  data. 
During  the  initial  testing  used  to  refine  the  techniques,  the  lenslet  size  was  480pm  and 
later  was  upgraded  to  a  higher  fidelity  MLM,  with  lenslets  of  133pm  across.  The 
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Wavescope®  software  allowed  various  inputs  for  each  test.  The  initial  testing  with  the 
480pm  MLM  was  completed  with  5  measurements  per  test,  5  frames  for  each 
measurement  collected  at  30Hz,  while  collecting  one  data  point  for  each  lenslet.  The 
finer  133  pm  MLM  produced  a  much  larger  quantity  of  data;  therefore,  the  tests  were 
reduced  to  3  measurements,  five  frames  per  measurement  at  5  Hz.  All  the  test  data  were 
exported  from  the  Wavescope®  software  to  data  files.  MATLAB®  script  was  written  to 
process  this  data  and  output  raw  displacement  data,  Zemike  plots  and  surface  smoothed 
plots. 

Mirror  Deflection  Results 

The  objective  behind  the  above  experiment  was  to  show  the  controllability  of  the 
polymer  membrane  through  voltage  differential  across  specific  piezoelectric  designs. 
There  was  definite  proof  that  the  0.25”  thick  membrane  and  piezoelectric  backing 
displaced  together  as  seen  in  the  non- actuated  and  actuated  displacement  plots  in  Figure  9. 


M2  TO  OVM  KMK  XOMI 


Figure  9.  Non-Activated  and  Activate  Membrane  Plots 
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For  the  most  part,  the  piezoelectric  placement  was  even  and  symmetric.  This  may  not  be 
the  case  for  a  final  membrane,  as  there  may  be  a  desire  to  activate  zones  of  different 
shapes  across  the  whole  membrane.  This  study  was  not  concerned  with  that  aspect,  but 
only  to  show  that  it  is  possible  to  displace  the  membrane. 

The  results  of  all  these  tests  supplied  plots  of  raw  data  of  the  surface  of  each 
membrane.  The  polymer  membranes  were  relatively  optically  flat  (3.95  X  PV,  0.63  X 
RMS);  however,  there  was  definite  high  frequency  distortion,  at  submicron  resolution, 
across  the  whole  membrane.  This  was  partly  resolved  by  “zeroing”  out  the  initial  0V 
deflection  plot  from  each  raw  data  plot  for  each  voltage  condition.  The  “zeroing”  process 
was  accomplished  by  subtracting  out  the  initial  wavefront  gathered  when  no  activation 
was  applied.  This  self- referencing  technique  leads  to  a  great  surface  flatness 
improvement.  For  example,  the  zero  voltage  case  when  self-referenced  is  nearly  an  order 
of  magnitude  better  than  the  raw  data  (0.59  X  PV,  0.07  X  RMS). 

Much  of  the  distortion  encountered  for  each  voltage  case  was  reduced  using  this 
method,  but  another  approximation  used  within  the  optics  community,  namely  Zemike 
polynomials,  supplies  another  way  to  look  at  the  data.  The  Zemike  polynomials  are  a 
way  to  express  the  wavefront  deviation  of  each  term  with  respect  to  the  actual  wavefront. 
In  other  words,  a  combination  of  a  series  of  multipliers  (Zemike  coefficient)  and  polar- 
represented  polynomials  (Zemike  polynomials)  tallied  to  give  a  representation  of 
displacements. 
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As  stated  before,  the  two  PVDF  membrane  configurations  have  the  important  test 
results.  There  were  numerous  tests  completed  on  both  membranes  at  different  voltage 
conditions.  Ml  was  used  mostly  to  determine  testing  procedures.  The  first  three  tests 
were  completed  before  the  second  layer  of  RTV615  was  applied.  The  following  figures 
are  deflection  plots  of  Ml  and  their  corresponding  Zernike  plots. 

After  applying  positive  or  negative  voltages,  it  was  noticed  that  there  was  some 
residual  deflection  in  the  region  of  actuation.  This  may  be  a  residual  charge  on  the 
control  surface.  The  WaveScope®  was  recalibrated  to  compensate  for  this  residual 
before  more  tests  were  conducted.  Much  more  was  realized  throughout  the  testing  of  Ml, 
but  the  most  important  aspect  was  that  these  tests  proved  low-order,  global  shape  control 
of  lightweight  membrane  mirrors  is  achievable  using  piezoelectric  control  layers,  such  as 
PVDF.  The  order  of  magnitude  of  this  control  is  the  same  as  that  of  the  surface  features, 
demonstrating  that  it  might  be  possible  to  correct  for  them. 

After  the  knowledge  gained  from  the  testing  of  Ml,  M2,  which  had  different 
control  zones  (refer  back  to  Figure  6)  was  tested.  As  with  Ml,  M2  was  evaluated  before 
any  voltage  was  applied  to  check  the  surface  flatness.  The  control  region  7  (a  3- mm 
circle  in  the  center)  tests  are  the  most  important  to  look  at  for  this  study,  since  they 
represent  an  axisymmetric  case.  Other  results  can  be  seen  here,  but  are  only  to  show  the 
ability  of  activation  and  displacement  of  other  areas  of  the  membrane. 

Control  region  7  was  first  activated  with  300,  followed  by  600V.  The  magnitude 
of  the  depression  increased  by  approximately  one  wave  as  the  voltage  was  increased 
(2.33  A,  versus  3.22  A  ,  respectively).  Again,  the  following  plots  are  the  displacement 
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plots  along  with  their  respective  Zernike  plots  for  each  case.  Figure  10  is  the  300  V  case 


and  Figure  1 1  is  the  600  V  case. 


M3-TV;  30CV  Sroctf-d  PV*  Itf-TIc  XBV /•tnfcn  C*4 


Figure  10.  300V  Experimental  Plots 


M7-TVS  env  VuttM  Art  V?-T  Id  Stv  fmrtkm  PM 


Figure  11.  600V  Experimental  Plots 


After  the  positive  voltage  tests,  the  control  leads  were  grounded  to  discharge  the 
PVDF  layer  and  retested  at  0V.  Figure  12  corresponds  to  the  0V  case.  Notice  just  a  little 
distortion  in  these  plots  caused  by  the  deflection  of  the  residual  positive  voltage. 
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M7-T 1*  C*  Irmaitmd  PM  W-TI#  IT,'  7mitt  Pw 


Figure  12.  OV  after  positive  voltage 


The  tests  were  re- accomplished  at  the  same  magnitude  negative  voltages.  There 
was  again  about  a  one-wave  change  in  the  deflection  (2.39  X  PV  for  -300V  and  3.4  X  PV 
for  -600V).  The  width  of  the  deformed  area  was  slightly  larger  than  the  control  region. 
The  deformed  area  was  about  40mm  and  50mm  for  -300v  and  -600V,  respectively. 

Figure  13  is  the  -300V  case  and  Figure  14  is  -600  in  the  center  region.  Notice  that  the 
negative  voltage  produced  a  positive  displacement.  The  Zemike  plots  seem  to  still  have  a 
representation  of  the  raw  data. 


VCM  n  -xm  r.o n* 


M7-TU  -2TCW  /»rn  kj>  Ftr 


Figure  13.  -300V  Experimental  Plots 
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M7-T«o  MM  lo  /wntu  PH 


Figure  14.  -600V  Experimental  Plots 

Once  more  the  0V  case  is  included  after  the  negative  voltages  are  applied  to  check 
for  any  effects  of  the  deformation.  Figure  15  show  that  there  is  a  trace  of  a  positive 
deformation,  caused  by  the  negative  voltage. 


M7>T  ih  CV  Smottal  P%*  UMIh  IN  7mrUm  Pw 


Figure  15.  0V  after  negative  voltage 

As  stated  before,  various  zones  were  activated  alone  and  together  to  verify  the 
controllability  of  the  membrane  surface.  They  are  included  here  only  to  show  the  effects 
of  these  zones  and  are  not  used  later  in  the  study.  Figure  16  shows  zones  3  and  6 
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activated  with  300  volts  in  each  section.  Notice  how  much  control  over  the  membrane  is 


possible.  Many  different  zones  could  be  designed  that  would  allow  a  fair  amount  of  local 
control  without  overly  affecting  the  adjacent  zones.  This  concludes  the  discussion  of 
work  completed  by  Sobers,  which  demonstrates  a  very  crucial  ability  of  static  control 
over  a  membrane  at  the  minimum. 


M7-T7t>  MW  Snathae  MT-fTb  XKV /am *»  Pfc* 


Figure  16.  300V  activation  in  Region  3  and  6 


Material  Property  Test  and  Results 

The  material  properties  of  many  polymers  are  unknown  due  to  some  complications  in  the 
testing  procedure  and  the  lack  of  need  in  the  past;  however,  the  Modulus  of  Elasticity  and 
Poisson’s  Ratio  are  required  components  to  the  finite  element  solution.  A  method  was 
therefore  designed  by  the  Air  Force  Research  Laboratory  and  some  affiliates  using  a 
combination  of  Moire  Interferometry  and  tensile  testing.  The  method  was  originally 
designed  for  testing  of  3-dimensional  woven  materials,  but  was  easily  adapted  to  the 
testing  of  the  RTV615  polymer  used  in  the  preceding  mirror  tests. 
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In  most  tensile  tests,  where  a  material  similar  to  a  metal  is  being  tested,  the  grips 
of  the  tensile  testing  machine  can  be  directly  attached  to  the  specimen  in  a  region  that  is 
not  to  be  involved  in  the  test.  This  can  be  accomplished  if  the  material  will  not  compress 
significantly  under  the  grip  strength  and  will  not  tear  as  loaded;  however,  this  would 
happen  with  a  polymer.  A  carbon  fiber  two-end  mounting  tab  system  was  designed  that 
would  serve  as  the  “grip  zones”.  They  actually  were  center  drilled  and  a  string  was 
inserted  to  hold  the  specimen  from  the  top  and  a  “bucket”  where  the  weight  would  be 
applied  on  the  bottom. 

The  preparation  of  the  polymer  specimen  involved  the  machining  of  a  few 
aluminum  pieces  about  2”  long  x  .25”  thick  to  be  used  as  a  mold  while  the  two  part  liquid 
polymer  RTV  615  was  poured  and  cured.  Aluminum  was  used  to  ensure  the  cured 
polymer  could  be  released  from  the  sidepieces  without  destroying  the  specimen,  since  the 
polymer  would  not  adhere  to  this  metal.  The  carbon  fiber  pieces  were  also  machined  to 
about  Vi”  x  ‘A”.  Two  aluminum  end  tabs  were  machined  to  about  the  same  dimensions 
and  a  transverse  hole  was  drilled  in  them  to  be  used  in  conjunction  with  thread  as  the 
mechanism  to  hold  the  ends  of  the  specimen.  The  specimen  can  be  seen  in  Figure  17. 
Notice  that  one  of  the  aluminum  end  tabs  has  fallen  off.  The  carbon  fiber  was  chosen  for 
its  porous  quality  into  which  the  polymer  could  infuse,  so  the  pieces  would  not  separate 
during  the  test.  After  accomplishing  the  machining,  a  diffraction  grating  was  selected 
that  would  serve  as  the  mechanism  to  imprint  onto  the  specimen  an  equivalent  tool  most 
closely  related  to  a  strain  gauge  in  a  normal  tensile  test.  This  diffraction  grating 
imprinted  polymer,  the  A1  end  tabs  and  carbon  pieces  are  the  only  parts  of  the  specimen 
tests.  All  other  parts  were  removed. 
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Figure  17.  RTV  615  Specimen 


The  specimen  was  created  in  the  following  fashion:  The  two  aluminum  and  two 
carbon  fiber  pieces  were  glued  onto  an  area  of  the  diffraction  grating  such  that  a  channel 
of  about  1.5”x  Vi'  was  created,  so  the  RTV  would  not  be  able  to  escape  from  the  channel. 
The  polymer  was  mixed  thoroughly  by  hand  and  placed  in  a  centrifuge  machine  to 
remove  all  the  air  pockets  and  bubbles.  During  this  time,  a  bit  of  polymer  bonding  agent 
was  applied  to  the  carbon  pieces  to  further  enhance  the  bonding  of  the  polymer  and 
carbon.  The  polymer  was  then  poured  into  the  channel  onto  the  diffraction  grating.  It 
was  allowed  to  cure  for  about  a  week,  removed  with  an  imprint  from  the  diffraction 
grating.  The  aluminum  pieces  were  also  removed  from  the  sides  of  the  specimen. 

The  next  step  was  securing  the  specimen  in  the  test  setup  by  threading  the  two 
holes  in  the  carbon.  The  specimen  was  hung  by  this  thread  at  the  top.  A  small  bucket  was 
hung  from  the  other  side,  toward  the  ground,  of  the  specimen  where  the  weight  would  be 
applied.  Different  attempts  at  constraining  the  loading  angle  were  made  to  ensure  little, 
or  ideally,  no  bending  when  weight  is  added.  Different  tests  were  accomplished  and  the 
early  tests  were  unsuccessful  mostly  due  to  bending  in  the  specimen.  New  constraints  of 
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the  loading  angle  using  bars  holding  the  thread  to  ensure  vertical  displacement  were 


applied  and  the  test  was  rerun. 

As  with  most  tensile  tests,  the  specimen  was  loaded  in  the  axial  direction  with  the 
procedure  stated  above.  Prior  to  loading,  the  diffraction  grating  was  brought  into  view 
through  the  Moire  Interferometer,  as  seen  below  in  Figure  18. 


Figure  18.  Moire  Interferometer 


This  is  accomplished  by  exposing  the  diffraction  grating  to  two  collimated,  coherent 
beams  of  laser  light.  A  schematic  is  shown  here  in  Figure  19.  Before  deformation  of  the 
specimen  is  accomplished,  one  diffraction  order  from  each  beam  is  assured  to  be 
perpendicular  to  the  specimen,  by  adjusting  the  angle  of  incidence  of  the  two  beams. 
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Figure  19.  Moire  Schematic 


The  specimen  is  deformed  causing  the  diffraction  grating  frequency  to  change  and 
warpage  appears  in  the  flat  wavefronts  diffracting  from  the  specimen.  An  interference 
fringe  pattern,  or  a  contour  map  of  in-plane  displacements  in  direction  perpendicular  to 
grating  lines,  is  the  result.  Another  set  of  incident  laser  beams  are  rotated  90  degrees 
about  the  z-axis  from  the  above  laser  light.  The  beams  in  the  y-z  plane  are  blocked.  This 
is  accomplished  to  measure  displacements  in  two  orthogonal  directions. 

The  displacements  are  obtained  through  a  simple  relation  as  seen  in  the  following 
equations: 

U(x,y)  =  jNx(x,y) 

V(x,  y)  =  ^-Ny(x,  y) 

f  (3) 

(Mollenhauer) 

where  f  is  twice  the  diffraction  grating  and  Nx  and  Ny  are  the  “fringe  orders”,  a  numerical 
value  assignment  to  each  fringe.  Due  to  the  small  deflections,  the  strains  may  be 
determined  through  standard  small  deformation- strain  relationships  of  the  derivatives  of 


28 


the  above  displacement  equations.  It  is  possible  to  perform  the  analysis  either  manually 


or  using  advanced  digital  phase- shifting  techniques.  For  this  analysis,  about  80%  of  the 
width  in  the  center  of  the  specimen  was  used  for  the  calculations.  The  strains  were 
calculated  using  a  standard  numerical  differentiation  scheme .  The  detailed  information 
corresponding  to  the  tests  can  be  seen  in  Table  1.  The  two  tests  were  averaged  to  obtain 
147  psi  (approximately  1  Mpa)  and  0.497  for  the  Elastic  Modulus  and  Poisson’s  ratio, 
respectfully.  These  numbers  are  the  starting  point  for  the  finite  element  modeling 
analysis.  The  development  of  the  finite  element  code  is  covered  in  the  next  chapter  and 
the  results  covered  in  the  second  half  of  Chapter  4. 


Table  1  Tensile  Test  Details 


Test  # 

Loading 

Axial 

Strain 

Transverse 

Strain 

Shear  Strain 

Poisson's 

Ratio 

Modulus 

1 

10.528 

1686 

-789 

-112 

0.468 

143 

2A 

7.018 

1025 

-540 

-57 

0.527 

157 

2B 

14.037 

2206 

-1156 

-106 

0.527 

146 

Average  2 

10.5275 

1615.5 

-848 

-81.5 

0.527 

151.5 

Average 

10.5278 

1650.75 

-818.5 

-96.75 

0.4975 

147.25 
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Ill:  Theoretical  Development  and  FEM  Code 


Introduction 

Now  that  the  ability  to  have  some  active  and  specific  control  over  transverse 
membrane  deflection  has  been  shown,  a  way  to  model  this  behavior  would  prove  useful 
to  discover  the  many  possibilities  of  control  over  the  membrane  without  having  to 
accomplish  physical  tests,  especially  when  the  membranes  are  very  large.  The 
commercial  market  of  finite  element  modeling  was  investigated  to  attempt  to  accomplish 
the  solution  of  a  piezo -silicone  membrane  mirror.  Unfortunately,  none  of  the 
investigated  software  packages  met  all  the  requirements  for  this  analysis;  Some  of  these 
requirements  included,  but  not  limited  to,  micron  level  displacements,  a  piezoelectric 
element  that  suited  the  purposes  of  the  experiment,  and  modeling  the  polymer  to  behave 
in  the  manner  of  a  membrane.  Therefore,  after  extensive  theoretical  development  in 
beam  and  membrane  finite  element  modeling,  a  new  mathematical  methodology  based  on 
fundamental  perturbation  techniques  was  born,  and  given  the  name  Method  of  Integral 
Multiple  Scales  (MIMS).  (Rogers,  2001) 

This  technique  has  selectable  precision  when  applied  to  a  Langrangian 
represented  class  of  dynamic  systems.  As  a  proving  ground,  MIMS  was  first  applied  to  a 
simple  linear  beam  and  was  shown  to  produce  boundary  layer  results.  The  method  was 
applied  fully  to  a  finite  element  approach  and  the  accuracy  was  shown  to  be  three  orders 
of  magnitude  greater  than  standard  finite  element  formulation.  This  FE  methodology  was 
finally  applied  to  nonlinear  beam  and  axisymmetric  circular  membrane  to  compare  to 
analytical  solutions.  (Rogers)  At  the  time  the  dissertation  was  completed,  the  objective 
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was  to  setup  the  solutions  so  they  matched  the  analytical  solutions,  but  was  not 
necessarily  representative  of  the  experimental  solution. 

Beam  String 

A  beam- string  is  represented  by  an  inner  and  outer  region  to  describe  the 
deflection  The  outer  region  or  the  center  (i.e.  away  from  the  boundary  condition)  of  the 
beam- string  is  represented  by  string  equations,  while  the  inner  or  edge  region,  is 
represented  by  beam  equations.  An  energy-based  derivation  method  is  used  produced  the 
necessary  equations  used  in  this  analysis.  The  details  of  the  total  derivation  are  not 
included  here,  but  an  outline  of  the  method  is  shown  to  familiarize  the  reader  with  the 
notation  and  to  setup  an  example.  The  literature  has  the  complete  derivation  in  Chapter  3 
of  Rogers.  (Rogers,  2001)  A  strain-energy  representation,  the  potential  energy  of  an 
elastic  beam  is  defined  by: 

F=j|{£}r{o}-{£}r{a0}  +  {£0}r{o}JR 

+  -^K0w2x(0,  t)  +  —  Klw2x(L,  t)  (4) 

where  £  and  s  are  the  strain  and  stress  terms,  w,x  terms  are  the  slope  of  the 
displacements  at  the  boundaries  and  the  K  terms  are  the  spring  stiffness  at  the  membrane 
boundary. 

The  beam’s  kinetic  energy  can  be  expressed  as: 

T=  jip  {ul  +  w))dV  (5) 

V  ^ 
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where  p  is  the  density  and  the  u  and  w  terms  are  the  time  derivative  of  the  displacements. 


Making  the  assumption  of  small  deflection  with  respect  to  overall  dimensions,  the 
system’s  non- conservative  forces  from  a  pressure  force  (P)  and  displacements  (u  and  w) 
can  be  approximated  by: 

W  =  -Jp(w  +  m  xw-  w Xu)dx  (6) 


Analysis  of  a  laminate  requires  further  analysis,  seeing  as  each  layer  may  have 
independent  material  properties.  The  laminate  could  be  reduced  to  a  one -dimensional 
integro- differential  equation  when  the  through-the-thickness  strain  is  assumed  to  be 
constant  at  any  cross  section.  Using  Hamilton’s  Principle, 


J  (ST  -SV  +  8W)dt  =  0 


(7) 


where  d  is  used  to  represent  variation.  The  following  dimensionless  variables  are  defined 
as: 


’  L 
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£4 
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EA  *  EA 
-q  EA,...  =  — 

1  Ez  LEA 


LEA 


(8) 


f2  _  C2  t2 


The  variation  is  carried  out  and  the  above  dimensionless  variables  introduced  to  produce 
the  following  equations  of  motion: 
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u,„  -TKu  =  -2om^  +^-(n  - 1  +  EA )^(vva  - 2m-X) 
- 2/3  vvi  -  ^  [' w,«  (!  -  | w5  +  M,xwi  ] 


(9) 


w,«  -0-£4  )w„  +  r2w  u.xx  =  -2(1  +  (n  - 1  +  £4  )— (ewx) 

ox 

+  P(\-  2u  x )  -  £A:  [m  q.  (1  - 1  w2  +  M  XW2  ] 


(10) 


e  = 


1 

+  — wr 

2  ’* 


(11) 


These  equations  must  satisfy  the  following  non-dimensional  boundary  conditions: 
u,  (0)  =  ut£ (0)  =  0  w/0)  =  0  N20w2%  (0)  =  -KQwn  (0) 

N20wm  (0)  =  - K0w 3^  (0)  -  EAu3  N20w4%  (0)  =  -K0w4A (0)  +  EAzlih:-  (12) 

where  i  =  3,  4,  5  and  j  =  2,  3,  4. 

At  this  point,  a  perturbation  solution  is  sought  on  a  radius  of  gyration,  r,  and  using 
the  following  non-dimensional  expansions : 

X=C0f  Tn  =  £  "X  03  =G30 +£03,  +£2G32  +... 

w(jC,X,£)  =£  2W2(x,X  )  +  £3W3(i:,X)  +  ... 

u(x,x;  £)  =e3u3(x,t  ) +e  4m4(i,x  )  + ...  (13) 

,  d 

D  =  Dn  +eD,  +£  D0  +...  £  =  r  D  = - 

o  1  2  ,  a2. 

The  string  equation  is  produced  by  substituting  the  expansion  into  the  differential 
equation  and  collecting  powers  of  £  .  The  coefficient  of  the  lowest  power  of  £  produced 
the  following  equation: 
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d~w2  1  —  EAe  jd“w2 


.  2iL  ^H+A 

C00  dx  COq 


After  solving  this  equation  for  W2,  it  is  noted  that  only  two  boundary  conditions 
are  satisfied.  This  is  not  enough  to  solve  the  unperturbed  equation  To  satisfy  the 
complete  set  of  boundary  conditions,  it  is  necessary  to  examine  the  original  equations  in 
the  neighborhood  of  the  boundary.  This  is  accomplished  by  rescaling  x  near  the 
boundary.  For  example,  E,  -  x/  e  is  how  this  rescaling  is  accomplished  at  x  =  0  and 
E,  =  ( L  -  x)/e  at  x=  1.  Notice  that  c,  is  zero  at  x=0,  while,  if  x  is  any  other  value,  c, 
approaches  infinity  as  £  approaches  zero.  Using  Equations  9-11,  and  following  the 
above  procedure,  the  resultant  is: 


3'w,  d4w , 


-  2  dw  (.  r.  \1  3  ..  2  1  x£2 

-2(ie  t — +(l-£AE)— — ((mx -Mf - — )  (15) 

oX0  £  dc,  2  £  oc, 


To  complete  the  solution,  the  behavior  near  the  boundary  must  be  matched  with 
the  membrane  solution,  using  Matched  Asymptotics  (Rogers,  2003).  This  identifies  the 
thickness  of  the  boundary  layers.  The  result  of  this  matching  leads  to  the  deflection 
results: 


uc{x)  =  0 

wc(x)=e2^-(x- x2)- £3[ro(i  —  *-  e~xl£)  +  Tx(x- e~~ )] 


where 


->  _  go  r^2  + 

°  K0-Nl  2  Kq 


r  -  r El  4-  ^.3] 

1  K,-N20  2  Kx 
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The  K  terms  are  the  spring  constants  defined  (with  ‘hats’)  in  Equation  8  and  characterize 
the  torsional  spring  boundary  conditions  at  the  ends  of  the  membrane.  No  is  the  tension 
also  defined  in  Equation  8.  EAe  Z3  represents  the  non-dimensional  moment  imposed  by 
the  piezoelectric  layers.  Now  that  the  static  solution  is  fully  developed  for  a  beam  string, 
an  example  solution  to  observe  the  effects  of  some  of  the  variables  (material  properties  or 
system  properties)  is  completed  here  before  continuing  onto  the  membrane  problem. 

This  is  a  mini-parametric  study  of  a  beam  string  problem  changing  a  few  of  the 
variables  to  see  the  effects.  These  numbers  will  be  similar  to  the  numbers  that  will  be 
used  in  the  parametric  study  of  the  membrane  solution.  The  code  for  this  problem  was 
written  by  Rogers  in  MathCAD.  The  importance  of  the  results  lies  in  the  magnitude  of 
the  deflection;  therefore,  the  plots  will  be  quickly  talked  about  for  each  case.  The  model 
is  represented  by  a  two  layer  system  (1  membrane,  1  piezo  layer);  however,  this  is  just  a 
demonstration  of  the  effects  of  material  properties. 

The  first  variable  to  be  varied  is  the  tension  term,  defined  by  No  in  the  equations. 
A  few  tension  values  tested  are  1,  5,  10,  100,  and  1000  N.  The  voltage  condition  used  is  - 
300  V  on  the  piezoelectric  layers,  no  pressure  term  thicknesses  of  0.0001  m  for  both 
layers,  2.8  GPa  and  1.8  GPa  for  the  membrane  and  PVDF  moduli,  respectfully.  These 
cases  can  be  seen  below  in  Figure  20- Figure  24,  in  the  order  defined  above.  Both  axes 
are  scaled  on  all  the  following  figures.  The  abscissa  is  scaled  by  the  length  of  the  beam 
(x=x/F)  and  the  ordinate  contains  the  radius  of  gyration  and  an  actual  displacement. 
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Figure  20.  Beam  String  Deflection  for  No=l 


Figure  21.  Beam  String  Deflection  (No=5) 


Scaled  Radial  Distance 

Figure  22.  Beam  String  Deflection  (No=10) 
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Beam  String  Deflection  Representation 


Scaled  Radial  Distance 

Figure  23.  Beam  String  Deflection  (No=100) 


Figure  24.  Beam  String  Deflection  (No=1000) 


Notice  how  the  increase  in  the  tension  term  decreases  the  deflection  magnitude. 
This  is  an  obvious  effect  of  the  tension  applied  at  the  two  boundaries.  The  other  effect  is 
that  the  beam- string  initially  behaves  as  a  string,  where  the  piezo  layer  is  deflecting  the 
structure  in  a  parabolic  shape,  with  the  boundaries  only  holding  to  the  fixed  condition. 

As  the  tension  is  increased,  the  structure  starts  to  become  stiffer  and  therefore  has  a  very 
small  deflection  that  occurs  across  the  whole  boundary,  with  a  greater  sbpe  at  the 
boundary. 
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The  next  parameter  that  is  investigated  in  this  study  is  the  thickness.  The  tension 
value  used  was  25  N  for  this  case.  It  does  not  have  any  specific  correlation,  but  was 
chosen  just  to  have  reasonable  results  for  a  few  of  the  thickness  conditions.  The  moduli 
were  kept  constant  as  defined  above  for  the  tension  study.  The  different  thicknesses 
tested  were  as  follows,  along  with  their  percentage  related  to  the  “original”  thickness: 
0.00001  m(l/10th),  0.0002  m(2x),  0.001  m(10x),  0.005  m(50x),  and  0.00635  m(.25”). 
This  study  was  based  solely  on  trial  and  error  and  therefore  these  thicknesses  were 
chosen  as  a  basis  for  what  might  hold  in  some  real  cases.  Figure  25 -Figure  29  show 
these  results. 


Figure  25.  Beam  String  Deflection  (t=le-5  m) 
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Figure  26.  Beam  String  Deflection  (t=2e-4  m) 


Beam  String  Deflection  Representation 


Figure  27.  Beam  String  Deflection  (t=0.001  m) 
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Beam  String  Deflection  Representation 


Figure  28.  Beam  String  Deflection  (t=0.005  m) 


,  Beam  String  Deflection  Representation 
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Figure  29.  Beam  String  Deflection  (t=0. 00635  m) 

Notice  the  quick  changes  as  the  thickness  increases.  This  study  shows  how 
important  it  is  to  know  the  thickness  and  how  it  effects  the  small  epsilon  term.  The  larger 
thickness  causes  the  epsilon  terms  to  grow  to  the  point  that  the  solution  can  no  longer  be 
considered  valid.  The  first  plot  for  the  small  thickness  possesses  very  small  displacement 
and  the  tension  seems  to  have  dominance.  It  behaves  as  a  membrane  with  the  small 
thickness  and  a  beam  as  the  thickness  is  increased  by  an  order  of  magnitude.  As  the 
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thickness  is  further  increased,  the  solution  starts  to  blow  up.  The  one  with  the  same  order 
of  thickness  has  a  reasonable  string- like  solution.  As  soon  as  the  thickness  becomes  too 
large,  the  result  stops  resembling  a  beam  string  and  has  no  real  meaning.  The  reason 
behind  the  negative  displacement  in  the  thicker  membrane  is  caused  by  the  solution 
blowing  up.  The  thicker  layer  seems  to  cause  the  solution  to  produce  results  of  large 
deflection  and  the  layer  pulls  the  whole  system  flat.  This  shows  the  sensitivity  to 
thickness  changes  in  the  non-piezo  layer  which  violated  the  small  perturbation  term 
defined  by  the  radius  of  gyration  The  thickness  seems  to  have  to  be  on  an  order  of  10‘4 
for  the  other  conditions  to  have  a  practical  result.  Once  the  thickness  was  increased 
enough,  the  solution  doesn’t  maintain  an  actual  result.  Now  that  the  system  properties 
have  been  evaluated,  the  next  step  is  to  test  the  effects  of  a  material  property. 

The  main  material  property  that  should  have  an  effect  on  the  deflection  is  the 
modulus  of  the  polymer  layer.  A  range  of  magnitudes  were  tested,  namely  106,  107,  108, 
and  109  Pa.  The  piezo  layers  were  maintained  at  1.8  GPa.  They  can  be  seen  below  in 
Figure  30- Figure  33. 
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Figure  30.  Beam  String  deflection  (E=106) 


Figure  31.  Beam  String  deflection  (F=107) 


Figure  32.  Beam  String  deflection  (E=108) 
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Figure  33.  Beam  String  deflection  (E=109) 


Somewhat  surprising  is  the  lack  of  sensitivity  as  the  modulus  is  increased.  This 
may  just  be  due  to  the  piezo  layer  dominating  the  solution.  More  investigation  will  be 
placed  on  changing  the  piezo  properties  in  the  membrane  solution;  however,  since  the 
modulus  did  not  affect  the  deflection  solution  by  much,  a  larger  range  of  moduli  will  be 
studied  to  see  what  would  happen  with  different  materials. 

While  all  of  the  above  material  properties  show  changes  in  the  deflection  order 
and  shape,  some  of  them  actually  affect  the  solution  in  a  different  manner.  Recall,  e  is  a 
small  non-dimensional  term  that  plays  a  crucial  role  in  the  perturbation  techniques.  This 
term  actually  appears  in  the  shape  functions  used.  A  large  £  changes  the  shape  functions; 
thereby,  causing  the  solution  to  become  unrealistic.  The  finite  element  solution  is  only 
valid  for  the  small  epsilon.  This  is  evident  as  the  thickness  was  increased  and  the  plots 
started  to  have  results  that  are  not  actually  representative  of  a  beam  string  solution. 
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This  is  a  short  introduction  on  the  beam  problem  and  proves  useful  in  looking  at 
the  next  step,  the  axisymmetric  membrane.  The  whole  derivation  can  be  investigated  in 
the  literature,  although  the  important  points  are  shown  here  (Rogers,  2001).  Table  2 
shows  all  the  different  tests  in  the  above  parametric  study 


Table  2  Parametric  Study 


Test 

name 

Membrane 

Modulus 

(GPa) 

PVDF 

Modulus 

(GPa) 

Membrane 

Thickness 

(m) 

PVDF  Thickness 
(m)-Both  layers 

Tension 

Voltage 

Tension 

2.8 

1.8 

0.0001 

0.0001 

1 

-300 

2.8 

1.8 

0.0001 

0.0001 

5 

-300 

2.8 

1.8 

0.0001 

0.0001 

10 

-300 

2.8 

1.8 

0.0001 

0.0001 

100 

-300 

2.8 

1.8 

0.0001 

0.0001 

1000 

-300 

Thickness 

2.8 

1.8 

0.00001 

0.0001 

25 

-300 

2.8 

1.8 

0.0002 

0.0001 

25 

-300 

2.8 

1.8 

0.001 

0.0001 

25 

-300 

2.8 

1.8 

0.005 

0.0001 

25 

-300 

2.8 

1.8 

0.0635 

0.0001 

25 

-300 

Modulus 

0.001 

1.8 

0.0001 

0.0001 

25 

-300 

0.01 

1.8 

0.0001 

0.0001 

25 

-300 

0.1 

1.8 

0.0001 

0.0001 

25 

-300 

1 

1.8 

0.0001 

0.0001 

25 

-300 

Plate  Membrane  Equivalent 

Now  that  the  beam- string  problem  has  been  developed  and  a  parametric  study 
completed,  the  next  step  is  to  show  the  development  of  the  model  of  a  thin  membrane¬ 
like  system.  This  can  be  accomplished  using  a  plate- membrane,  where  low  bending 
stiffness  results  in  localized  effects.  Again,  a  Langrangian  system  is  used  for  this 
development.  The  equivalent  equations  related  to  the  beam- string  problem  will  be 
defined  here.  (Rogers,  2001) 
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The  potential  energy  is  developed  in  a  cylindrical  coordinate  system  of  r-0  -z  and 


is  seen  defined  here: 

F=ji{£f{o}-{£f{O0}+{£0}r{o}dl/  (18) 

v2 

+  jK0(6)w2r(R,t) 

The  K  term  represents  a  torsional  edge  spring,  which  was  incorporated  into  the  equations 
for  need  if  there  is  a  flexible  boundary.  The  potential  energy  can  be  rewritten  in  terms  of 
the  strain  relations  for  the  laminate.  The  kinetic  energy  can  be  represented  by: 


J^P  (u)  +  v]  +  w))dV 


(19) 


Assuming  constant  through- the-thickness  strain,  the  system  can  be  collapsed  into  a  one 
dimensional  integro- differential  system.  The  resulting  equations  of  motion  become: 


P hu„  —  EH —{[r(e r  +t)<?e)J  r  +  -J- (1  -1) )erQ e  -  ee  -Ver}--[Nr- Ng]  = 
r  2  r 

EH  - wm  +-wr  -  [ rw  +X>  - wee  ],  “  (1  -U )[- w,r9  -  \  w8  ],0  } 

y*  y  y  7*  7* 


(20) 
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(21) 


P  hv„ -EH-{[(ee  +\)e>,)]0  +^-(1 -'l))[e(0  +(re#)r}  = 
r  2 

EHz-{[\wm  +-w  +Dw  Je  +(l-D)[-wr0  ~\wp  +(-wr9  — l wB )  ] } 

y  y  y  y  y  y  y 

P hwtt  —Nr  —(rw  )  -  Nd  \ wee  -  2 N rQ -wrQ  +  DV4w  = 
r  r  r 

EH-{[rwr(er+X)ee)]r  +[-w0(e0  +  t)e,.)]0  +^(l-'U)[(w0e,.0 )  r  +  (w,.e,0 )  0  ]} 
r  r  2 

£Z-{[(e0  +ver)  +  rw r[wrr  +'D(-wr +4-wee )] +  (l-l))-we (w^,  --we)]r  (22) 

y  1  1  y  1  y  5  y  5  ’  y  ’ 

+  [(l-D)-erH  +-We[-wr  +\wm  +\)wrr]-(l-'D)-wr(wr0  --w0)]0 

y  y  y  y  y  y 

-[r(er  +X>ed)]rr  +  [-(ed  +Der)]ee  +(l-'uK0r0} 
r 

The  following  boundary  conditions  are  applied  to  obtain  a  solution: 

u(R,Q)  =  0  u  r(R,Q)  =  u  e(R,Q)  =  0  u(r,Q)  =  0 

v(R,  9)  =  0  v  r(R,Q )  =  ve  (R,Q )  =  0  v(r,9)  =  v(r,0  +2n) 

vj9(r,0)  =  v0(r,e  +2k) 

w(0,9 )  =  a;a  <  °°  wr(O,0)  =-wr(O,0  +7t )  (23) 

D[rw  JR,Q)  +D(- wee  (R,0 )  +  wr(i?,0))]  = 
r 

Nrz  -  EZ£  (i?,0 )  +  K(d  )wr(R,Q ) 

+  £Zr[(M ,.(/?, 0 )  + ^ w; (i?,0 )) -t) (-n(i?, 0 )  + i v 0  (i?, 0 )  + w20  (i?, 0 )] 

2  ’  r  r  ’  2  r 

These  systems  of  equations  are  put  in  non-dimensional  form  before  applying  the 
expansions  just  like  the  beam  string: 

X  =G3f  Tn=e\  03  =G3() +8  03,  +£2G32  +... 
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w(r,0,x;£)  =  £3w3(r,0  ,x)  +  £  4w4(r,0,x )  + ... 
v(r,0  ,x  ;e)  =  e  3v3(r,0  ,x )  +  £  4v4(r,0  ,x )  + ... 
w(r,0  ,x;£)  =  £2w2(r,0,x)  +  £3w3(r,0,x)  +  ... 


(24) 


D  —  Dq  +  eD  |  +  £  D  2  + . . .  E  —  r 


1 

T1  =-Tli 
£ 


From  this  point,  axisymmetric  conditions  can  be  applied,  thus  removing  all  0 
dependence.  It  is  also  assumed  EHm  =  EHer  =  EH,  ,  meaning  bidirectional  piezoelectric 

materials  are  being  used,  which  is  also  true  for  the  PVDF.  The  system  is  represented  by 
membrane  equations  in  the  center  of  the  system,  while  near  the  edges  they  are  beam 
equations  (boundary  layer).  Once  all  this  procedure  is  completed,  the  time  dependency 
can  be  removed  to  yield  a  closed  form  static  solution,  much  as  above  in  the  beam  string 
problem: 


wc(r )  =£2  — (1  -  r2)  +£3|  — 

4  [A'-l 

+  £ 4  j-^-TV 2  (1  +\)  )EHe  2 ( 1  -  r2  )| 


fL+EZi 
2  K 


^-0 


+  £ 


K 


K- 1 


~2K-3  +  2v 

(P2+EZeA 

EZe2 

Ts-i| 

1 

R-l 

{2  K  J 

K 2  _ 

(25) 


The  £  term  for  this  case  needs  to  remain  very  small  (<.01  for  most  cases).  The 
static  displacement  for  a  thin  membrane  plate  is  now  fully  developed.  The  finite  element 
approach  was  also  developed  by  Rogers  (2001);  however,  it  will  not  be  discussed  here, 
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but  the  shape  functions  are  scaled  in  a  similar  fashion  to  the  previous  parameters  as  a 
function  of  £  .  The  Langrangian  expansion  is  also  grouped  into  an  £  -order.  Applying 
Euler’s  equations  to  the  selected  Langrangian  element  produces  the  system  of  finite 
element  equations  that  are  directed  to  the  desired  solution.  Rogers  and  Agnes  (2003) 
have  performed  a  complete  stiffness  matrices  development  applying  C1  shape  functions. 
The  £  function  is  quickly  noticed  to  be  a  dominating  player  in  the  finite  element  solution 
and  thereby  kept  very  small  for  accurate  results. 

More  specifics  about  the  program,  as  well  as  parts  of  the  code  are  included  in  the 
Appendices.  As  stated  before,  the  above  theoretical  development  is  a  synopsis  of  the 
work  completed  under  another  study  and  only  serves  as  a  guide  for  how  development  of 
the  equations  for  a  membrane  was  created  and  the  steps  needed  to  program  a  new  finite 
element  code  that  can  handle  the  small  displacements  involved  with  a  mirror  system  as 
defined. 
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IV :  Results  and  Discussion 


Introduction 

Now  that  the  development  of  the  analytical  solution  to  the  axisymmetric 
membrane  plate  has  been  shown  and  the  finite  element  method  defined  by  Rogers 
discussed,  it  is  time  to  run  solutions  for  the  membrane.  The  beam- string  example  run  in 
the  previous  chapter  illustrates  the  need  to  perform  a  parametric  study  based  on  some  of 
the  material  and  system  properties  throughout  different  ranges.  This  parametric  study  is 
unrelated  to  the  experimental  study  from  Chapter  2;  however,  the  study  that  follows  this 
compares  the  finite  element  solution  with  the  experimental  results.  It  is  designed  to 
demonstrate  how  each  material  or  system  property  affects  the  solution  independently. 

The  input  properties  (material  and  system)  may  affect  the  solution  drastically  and 
therefore  a  parametric  study  was  designed  to  investigate  these  effects. 

Another  study  completed  is  the  comparison  of  the  reduced  experimental  results 
from  Sobers  and  the  finite  element  solution  based  on  realistic  material  properties.  This 
study  is  to  look  at  the  capabilities  of  the  finite  element  program  in  comparison  to  actual 
real  world  results.  This  can  serve  as  an  avenue  to  find  the  strengths  and  limitations  of  the 
program  to  see  what  can  be  completed  in  the  future  to  improve  applicability.  Before 
proceeding  into  the  results  of  these  studies,  a  short  explanation  of  how  they  are  laid  out 
follows. 

Zernike  Coefficients 

As  stated  numerous  times  throughout  this  thesis,  the  deflection  of  the  membrane 
is  very  small.  In  the  optical  community,  it  is  common  to  describe  the  wavefront  shape, 


49 


the  deviation  from  the  planar  and  perpendicular,  in  terms  of  Zemike  polynomials,  a  polar 


coordinate  oriented  system.  Through  the  use  of  a  multiplier  and  a  polynomial  with  radial 
and  angular  description  of  each  location  on  a  wavefront,  the  wavefront  deformations  can 
be  represented.  These  polynomials  are  unique  and  possess  desirable  properties  because 
of  their  orthogonality.  This  provides  the  opportunity  to  sum  the  terms  in  the  following 
fashion  for  the  wavefront  (W): 

L 

W(p,0)  =  £^£/r(p,0)  (26) 

r 

where  A  is  the  Zernike  coefficient  and  Ur  represents  the  Zemike  polynomials. 

p  represents  the  distance  from  the  center  of  the  membrane,  while  0  is  the  angle 
from  the  positive  y-axis  moving  in  the  clockwise  direction.  Many  texts  have  defined  the 
Zemike  polynomials  in  different  orders.  Table  3  below  is  a  portion  of  the  35  Zemike 
numbers  used  in  the  experimental  thesis  described  in  Chapter  2.  The  Zemike 
polynomials  used  in  the  finite  element  code  in  Chapter  3  were  ordered  in  a  slightly 
different  way.  This  was  not  a  concern  since  it  is  the  combination  of  all  Zemike 
polynomials  that  defines  the  wavefront;  therefore,  when  the  same  polynomials  are 
included  in  the  final  formulation  the  result  is  the  same. 


Table  3 :  Zernike  Polynomials 


Zernike  number 

Zernike  Polynomial 

Meaning 

1 

P  COS0 

Tilt 

2 

p  sin0 

Y  axis  tilt 

3 

2  p  2-1 

Defocusing 

4 

p  2cos20 

Astigmatism  with  axis  at  +/-  45  deg 

5 

p  2sin20 

Astigmatism  with  axis  at  +/-  0  or  90  deg 

6 

(3  p  2-2)  p  cos0 

Primary  coma  along  y  axis 

7 

(3  p  2-2)  p  sin0 

Primary  coma  along  x  axis 
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8 

6  p  4-6  p  2+1 

Primary  spherical  aberration 

9 

p  3cos(30  ) 

Triangular  astigmatism,  base  on  y  axis 

10 

p  3sin(30 ) 

Triangular  astigmatism,  base  on  x  axis 

11 

(4  p  2-3)  p  2cos20 

Secondary  Astigmatism  x 

12 

(4  p  p  2-3)  p  2sin20 

Secondary  Astigmatism  y 

13 

(10  p  4-12  p  2+3)  p  p  cos0 

Secondary  coma  x 

14 

(10  p  4-12  p  2+3)  p  sin0 

Secondary  coma  y 

15 

20  p  D-30  p  4+12  p  2-1 

Secondary  Spherical 

Reduced  Experimental 

As  explained  in  Chapter  2,  the  WaveScope®  exports  raw  data  vectors  and 
calculates  the  Zernike  coefficient  for  each  polynomial  based  on  an  approximation 
program.  The  data  was  exported  to  MATLAB®,  which  was  used  to  create  the  figures. 
Figure  34  is  a  sample  of  the  Zernike  plots,  both  in  3-D  plot  and  in  a  contour  plot  for  a  300 
volt  differential  in  the  center  3mm  of  the  membrane. 


M7-Ttc  Xtt *  7«n«  tv> 


Figure  34.  300V  Zernike  Plots 

The  wavefront  produced  from  the  raw  data  can  be  seen  here  in  Figure  35.  Notice  the 
similarities  in  the  general  shape  the  raw  data  and  Zernike  plots. 
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M?-Ttc  300  Mk 


Figure  35.  300V  Raw  Data 


As  seen  in  the  surface  plots  of  the  Zemike  and  raw  data  above,  it  is  seen  that  the 
polynomial  fit  for  35  Zernike  numbers  doesn’t  cover  all  the  small  peaks  and  valleys  that 
is  seen  in  the  raw  data;  however,  the  general  shape  (center  valley  shape  and  size)  is 
captured,  which  is  enough  to  describe  the  wavefront  for  most  applications .  The  Zemike 
plot  could  therefore  be  considered  an  accurate  representation  of  the  mirror  wavefront. 
The  majority,  30  of  35  to  be  exact,  of  the  Zemike  polynomials  are  not  symmetric  terms, 
as  a  result  of  the  angle  dependence.  It  was  determined  early  on  in  the  experimental 
problem  to  define  the  piezoelectric  actuation  zones  in  a  symmetric  layout  with  other 
zones  that  could  be  activated  with  asymmetric  results.  This  was  seen  in  the  7  actuation 
zones  defined  in  Chapter  2  for  M2.  The  center  zone  7  was  activated  for  most  voltage 
conditions  and  axisymmetric  solutions  were  therefore  produced. 

Unfortunately,  the  solution  still  had  some  asymmetric  effects;  however,  the  terms 
corresponding  to  these  30  terms  as  a  whole  happen  to  be  a  smaller  order  than  the  5 
symmetric  terms.  This  can  be  seen  first  in  Figure  34  above,  as  most  of  the  concentration 
on  the  contour  plot  is  in  the  center  section,  the  zone  activated.  There  is  some  deflection 
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in  other  areas,  but  as  can  be  seen  in  the  reduced  Zernike  contour  plot  in  Figure  36,  the 
reduced  data  is  representative  of  the  full  Zernike  data  set.  Obviously,  the  two  contour 
plots  are  not  exactly  the  same,  implying  that  there  is  some  non-axisymmetric  solution  in 
the  above  case.  This  was  known  previously,  as  is  evident  in  the  Zernike  coefficients  of 
the  30  non-axisymmetric  terms;  however,  the  effects  of  the  30  terms  need  to  be 
investigated  to  see  how  large  they  are  in  relation  to  the  other  5  terms. 


Figure  36.  Symmetric  Contour  Zernike  Plot 


This  investigation  was  completed  using  the  finite  element  code  to  plot  out  a  few 
cases  of  the  asymmetric  solution.  This  plot  is  shown  here  in  Figure  37.  On  all  the 
following  plots  the  center  of  the  membrane  is  on  the  right  hand  side.  This  applies  from 
now  until  the  end  of  the  paper.  The  plot  corresponding  to  the  symmetric  solution  is 
included  to  show  the  difference  in  order.  Recall  the  ordinate  is  defined  in  wavelengths  of 
light  (value  multiplied  by  633  nm).  These  solutions  provide  an  insight  and  create  some 
questions  into  where  the  difference  in  the  wavefront  comes  about. 
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Figure  37.  Non- Axisymme trie  deflections 

Figure  37  represents  the  experimental  non- reduced  solution  for  all  35  Zernike 
numbers  and  the  30  non- symmetric  polynomials  for  a  number  of  radial  lines  at  different 
angles  from  the  y  axis.  Both  of  these  are  included  to  show  how  much  of  an  effect  the 
non- symmetric  polynomials  add  (bottom  plot)  in  comparison  to  the  whole  deflection  (top 
plot).  The  bottom  plot  has  approximately  half  the  order  of  the  top  one.  The  value  of  the 
30  asymmetric  Zernike  polynomials  is  a  fraction,  albeit  a  large  fraction,  of  the  symmetric 
solution,  so  it  shows  that  eliminating  these  terms  is  quite  a  leap  of  faith.  An  asymmetric 
solution  comparison  would  prove  to  be  more  useful,  but  the  capabilities  of  the  finite 
element  code  only  allow  axisymmetric  solutions;  therefore,  the  experimental  results  have 
to  be  reduced  to  this  level  of  axisymmetric  to  be  checked  against  the  finite  element 
solution. 

Figure  38  is  the  5  symmetric  Zernike  polynomials  plotted  out  to  see  the  order  of 
the  center  deflection  in  comparison  to  the  non- symmetric. 
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I 


Experimental  Reduced  ’Symmetric)  Solution 


i  II  : 
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Figure  38.  5  Symmetric  Deflections 

It  is  noticed  that  the  effect  of  these  terms  is  not  completely  insignificant,  causing 
large  fractions  of  a  wavelength  in  wavefront  difference.  Since  this  information  comes 
from  an  experiment,  there  could  be  a  multitude  of  factors  that  lead  to  this  wavefront 
shape.  One  idea  for  this  effect  may  lie  in  the  initial  wavefront  caused  by  manufacturing. 
The  casting  process  was  not  performed  under  the  most  ideal  conditions,  i.e.  not  in  clean 
room  conditions  or  the  process  of  removing  the  air  bubbles  with  a  dental  pick.  This 
effect  may  be  shown  in  a  0  V  case,  as  shown  in  Figure  39.  The  scale  for  this  zero  voltage 
case  is  very  small  (10~3  X ),  so  the  solution  is  nearly  zero  everywhere,  but  it  is  noticed 
that  it  isn’t  exactly  zero  at  the  edge.  Remember  that  these  plots  for  the  experimental 
correspond  more  to  the  wavefront  and  not  just  the  deflection,  so  this  part  of  the 
membrane  was  actually  a  hill  compared  to  a  flat  wavefront. 
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Experimental  Reduced  (Symmetric)  Solution 
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Figure  39.  0V  Reduced  Experimental  Deflection 

Sobers  had  taken  action  to  avoid  this  effect  by  “subtracting”  out  the  0  V 
information  from  each  of  the  plots.  Unfortunately,  it  was  shown  by  Sobers  that  there  are 
residual  effects  from  applying  voltages,  so  there  cannot  be  any  guarantee  that  all  the 
effects  of  the  material  is  eliminated  in  the  solution. 

Another  possibility  resides  in  the  layout  of  the  piezoelectric  zones.  Sobers 
mentioned  Ml  may  have  had  a  “leak”  from  one  activation  zone  to  another.  M2  didn’t 
seem  to  have  this  effect;  however,  the  leads  to  the  center  zone  and  other  zones  may  affect 
the  solution  to  create  non-axisymmetric  deflections.  The  other  component  of  this  idea  is 
that  the  membrane  is  a  continuous  surface  and  the  axial  effects  of  the  bi-directional 
piezoelectric  layer  may  actually  cause  some  deflection  that  changes  other  areas  of  the 
membrane  in  a  non-axisymmetric  way.  The  last  part  of  this  effect  is  the  process  of 
etching  the  piezoelectric  layer,  meaning  they  may  not  have  been  perfectly  etched; 
therefore,  causing  some  leak  over  or  shape  change  in  the  membrane. 
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While  these  above  ideas  may  not  be  the  only  reasons  or  may  not  be  the  actual 
reason,  they  demonstrate  that  when  dealing  with  very  small  deflections,  many  factors 
may  lead  to  errors  in  the  solution.  It  may  be  possible  and  necessary  to  ignore  some  of 
these  errors.  The  motivation  behind  carrying  out  this  exercise  was  to  show  that  the 
axisymmetric  polynomials  provided  a  representative  solution  to  the  full  solution.  With 
this  being  the  case  and  the  reduced  solution  produced,  the  next  step  of  employing  the 
finite  element  code  could  be  prepared  to  process  different  conditions  that  could  be  then 
compared  to  the  experimental;  however,  as  stated  before,  the  material  properties  are 
crucial  to  the  solution. 

Parametric  Study 

The  finite  element  code  has  been  shown  to  be  very  sensitive  to  the  £  term  and 
can  therefore  have  solutions  that  are  no  longer  representative  of  the  analytic  solution  or 
the  experimental.  This  parametric  study  will  investigate  how  the  solution  responds  to 
various  changes  in  the  input  properties,  namely  the  boundary  tension  (N),  thickness  (t), 
Modulus  of  Elasticity  (E)  and  element  size,  determined  by  the  number  of  elements  used. 
The  membrane  is  assumed  to  have  a  constant  radius  throughout  all  of  these  studies  of  3”, 
which  is  representative  of  the  experimental  membrane. 

The  first  variable  tested  is  the  tension  at  the  boundary.  This  term  is  a  non- 
dimensional  force  term  represented  by  pulling  on  the  membrane  on  the  outside  boundary 
and  in  reality  is  caused  by  the  stretching  of  the  membrane  described  in  Chapter  2.  This  is 
an  unknown,  unmeasured  value.  There  may  be  ways  to  obtain  this  value  through  a 
vibration  test  or  other  method,  but  for  sake  of  this  study  this  was  not  completed.  During 
this  test,  the  first  membrane  simulation  to  be  tested  is  the  case  for  a  Kapton-PVDF  layup, 


57 


which  is  not  the  membrane  that  was  experimentally  tested,  but  serves  as  a  benchmark  to 
test  various  properties.  The  Kapton  has  material  properties  of  406  ksi,  0.006”  and  0.3  for 
modulus,  thickness  and  Poisson’s  ratio,  respectfully.  The  PVDF  is  modeled  with  261  ksi, 
0.003”  and  0.3  for  the  same  properties.  Nine  elements  across  the  radius  were  used, 
giving  an  element  size  of  approximately  1/3”.  A  300  V  differential  is  supplied  to  an 
approximate  3  cm  center  “circle”  in  the  PVDF  layer. 

Interestingly,  the  tension  has  a  major  effect  on  the  total  solution.  A  small  tension 
of  0.1  produced  a  center  deflection  of  7000  X  or  roughly  5  mm,  which  would  be  a  visible 
deflection  and  not  reality.  Figure  40  shows  this  deflection  plot  for  this  condition  of  N. 
Observe  the  negative  deflection  represented  by  this  N  value.  This  is  one  of  the  quirks  of 
the  tension  that  will  be  discussed  shortly. 
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Intuitively,  on  the  other  end  of  the  scale  (a  very  large  N),  the  solution  should  and 
does  show  very  small  deflections.  This  is  demonstrated  for  an  N  equal  to  109.  For  all 
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Figure  40.  Membrane  Deflection  (N=0.1) 
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intents  and  purposes,  the  deflection  produced  by  this  N  is  nearly  0  (order  is  10"6  A, ).  For 
this  large  tension,  the  deflection  is  now  positive,  as  seen  in  Figure  41. 
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Figure  41.  Membrane  Deflection  (N=109) 

The  change  in  deflection  sign  may  be  explained  in  the  fact  that  the  tension  is  in  a 
few  of  the  non-dimensional  terms,  therefore  mathematically,  there  could  be  a  point  where 
the  tension  can  change  the  non-dimensional  terms  enough  to  cause  the  deflection  to  flip. 
Therefore,  the  next  part  of  this  study  is  investigating  this  phenomenon.  A  few  values 
were  chosen  to  have  both  a  negative  and  positive  deflection.  Just  a  few  representative 
cases  are  displayed  here.  The  following  Figure  42  represents  tensions  of  2500  on  the  top 
and  10,000  on  the  bottom 
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Figure  42.  Membrane  Deflection  (N=2500, 10,000) 

These  are  shown  mainly  because  they  have  the  same  magnitude  of  deflection,  but 
are  opposite  in  sign.  Therefore,  there  is  definitely  a  snap  through  point  in  between  these 
two  values. 

Before  investigating  the  actual  point  of  this  happening,  a  few  more  values  were 
run  to  see  if  another  “singularity”  existed  to  produce  another  result  of  snap  through. 
Although  not  every  case  was  investigated,  enough  runs  were  completed  to  sho  w  that  once 
above  the  initial  snap  through  point,  up  to  the  109  case,  the  deflection  never  again  became 
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negative.  This  may  be  seen  in  Figure  43  for  the  tension  equal  to  25,000  and  the 
deflection  starts  to  become  small  as  it  approaches  very  small  values  for  much  larger  N 
values.  This  phenomenon  is  purely  a  mathematical  one  that  comes  about  from  the  non- 
dimensionalzing  of  the  input  variables.  There  is  not  likely  to  be  a  physical  representation. 
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Figure  43.  Membrane  Deflection  (N=25000) 


The  snap  through  point  is  of  some  interest,  since  it  represents  some  sort  of 
singularity  in  the  code  that  might  warrant  further  investigation,  if  the  code  is  ever 
modified  or  enhanced.  A  manual  convergence  study  was  completed.  A  middle  point 
between  2500  and  10,000  was  chosen,  the  program  run  a  few  times  until  it  was  noticed 
that  the  deflections  were  becoming  very  small.  One  more  case  was  run  at  N  equal  to 
3242.5.  Figure  44  shows  this  deflection  plot. 
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Figure  44.  Membrane  Deflection  (N=3242.5) 

The  piezoelectric  element  in  the  program’s  code  and  the  non- dimensioning  of  the  terms 
result  in  the  flip  effect  and  this  dip  at  this  value.  The  tension  at  the  edge  should  reduce 
the  deflection  as  it  is  increased,  but  it  doesn’t  seem  like  it  should  cause  an  opposite 
deflection.  The  deflection  scale  on  this  plot  is  on  the  order  of  10‘5  A ;  therefore,  this  is 
essentially  no  deflection.  The  dip  negative  before  the  upward  deflection  is  probably  just 
an  anomaly  that  is  inherent  in  the  code.  An  intermediary  study  at  this  point  may  prove  to 
yield  some  information  on  the  tension  in  regard  to  some  of  the  other  material  properties. 

Would  a  change  in  the  thickness  of  the  PVDF  layer  or  the  modulus  of  either  layer 
change  the  point  where  this  deflection  flip  happens?  The  thickness  of  the  PVDF  layer 
was  reduced  to  3xl0"5  in  and  this  created  changed  the  point  of  flip.  This  order  of  the 
thickness  is  more  realistic,  as  will  be  seen  in  the  actual  comparison  to  experimental 
results  later  in  this  chapter.  Increasing  the  PVDF  modulus  to  the  order  of  109  has  a 
similar  effect  of  eliminating  the  anomaly.  This  goes  to  show  that  the  non-dimenional 
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terms  are  easily  modified  to  change  the  flip  value  of  tension.  Since  the  modulus  and 
thickness  have  some  control  over  the  solution,  they  warrant  further  study.  Both  plots  can 
be  seen  in  Figure  45. 
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Figure  45.  Anomaly  Plots  (t=3x!0's  and  K|*vm  =  1 09> 


This  unknown  quantity  of  tension  has  quite  a  bit  of  control  over  the  final  solution. 
With  that  said,  there  are  ranges  of  N  values  that  would  make  more  sense  than  others. 
Obviously,  the  small  value  near  zero  and  the  very  large  value  produce  poor  results  and  do 
not  have  physical  significance,  i.e.  nearly  no  tension  for  N  equal  to  0.1  and  a  tension 
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where  the  epoxy  would  fail  for  the  109  case.  The  middle  case  is  also  not  an  acceptable 
tension;  therefore,  tension  could  almost  be  considered  a  subjective  variable  that  could  be 
modified  slightly  to  produce  reasonable  results.  For  most  ranges  of  tension,  the  general 
shape  of  the  deflected  membrane  was  the  same;  consequently,  N  may  be  treated  as  if  it 
were  a  linear  scale  factor  in  many  cases.  Not  having  a  secure  value  for  N  is  not  an  ideal 
situation,  of  course,  but  there  may  be  a  way  to  obtain  a  reasonable  value  for  it  as  the 
results  are  compared  to  the  experimental  later  in  this  chapter.  Experiments  may  be 
designed  to  test  and  discover  this  number  in  the  future. 

The  tension  does  play  a  major  part  in  controlling  how  much  deflection  occurs  in 
the  membrane,  but  it  is  not  one  of  the  variables  that  controls  the  critical  factor  of  £ 
mentioned  before.  This  is  directly  influenced  by  the  thickness  of  the  membrane  layer. 
This  study  investigates  the  changes  in  that  thickness  within  a  relatively  tight  range.  All 
properties  are  kept  the  same  as  the  tension  study  above  and  a  few  tensions  were 
investigated  for  each  thickness.  The  thickness  tested  were  +/- 10,  30  and  50%  from  the 
original  value  of  0.006”  or  0.0066,  0.0078,  and  0.009”,  respectfully. 

The  solution  in  this  range  of  thicknesses  was  very  sensitive.  The  best  example  of 
this  was  for  a  tension  that  was  small  (N  =  100).  Figure  46  below  is  this  case  for  0.0066” 
and  0.009”.  As  it  is  seen,  the  smaller  thickness  produced  approximately  35  X  of 
deflection  and  the  thicker  membrane  has  over  5000  times  that  displacement.  This  has  to 
be  caused  by  the  perturbation  epsilon  value  growing  too  large  for  the  integral.  The 
effects  of  all  the  terms  together  lead  to  the  integral  solution  to  blowing  up. 
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Figure  46.  Membrane  Deflection  (t=0.0066  and  0.009) 

For  other  values  of  tension,  the  result  wasn’t  quite  that  much  difference,  but  there 
still  was  still  on  the  order  of  4-5x  for  N=1000  and  5000.  Besides  the  scale  in  this 
thickness  range,  there  isn’t  much  else  to  change;  therefore,  the  other  plots  at  different 
tensions  values  are  not  included  in  the  paper.  While  these  thicknesses  are  not 
representative  of  the  experimental  membrane,  they  are  within  one  order  of  magnitude  of 
the  thickness  of  other  membranes  that  are  being  applied  to  optics  project  in  other  areas  in 
the  US  AF  (Marker,  2002). 
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Since  the  £  term  is  directly  related  to  the  thickness,  it  makes  sense  the  solution 
remains  sensitive  to  the  changes  in  thickness,  since  the  perturbation  technique  using 
MIMS  diverges  with  increasing  £  .  This  is  the  extent  of  the  study  dealing  with  the 
thickness  until  later  when  it  will  be  compared  to  the  experimental  results. 

Up  to  this  point,  only  system  properties  have  been  investigated.  They  are 
quantities  that  could  be  changed  through  methods  other  than  changing  the  actual  material. 
What  effect  does  the  material  have  on  the  solution?  There  are  basically  two  material 
properties  that  can  be  modified  for  the  static  solution,  namely  Young’s  Modulus  and 
Poisson’s  ratio.  The  next  part  of  this  study  investigates  different  moduli,  representing 
changing  the  material  being  tested.  During  this  test,  the  material  may  not  be  related  to  an 
actual  material,  but  is  designed  as  method  to  see  the  effects  on  the  solution. 

The  tests  involved  here  basically  used  a  range  of  E  values  from  100  to  106  psi  for 
the  membrane.  The  modulus  of  the  PVDF  layer  was  given  three  different  values  for  each 
value  of  membrane  moduli.  All  other  material  and  system  properties  were  kept  the  same 
as  the  original  tension  tests  and  the  tension  was  given  a  value  of  1000.  The  first  test  was 
using  a  PVDF  layer  with  261  ksi  for  its  modulus.  This  test  produced  results  that  behaved 
similar  to  the  change  in  tension.  At  one  point,  the  solution  flips  from  negative  to  positive 
deflection  and  there  happens  to  be  a  solution  where  there  is  a  dip,  similar  to  the  case  for 
tension.  The  solution  flips  near  an  E  value  of  240  ksi,  with  the  dip  noticed  at  239  ksi. 
Figure  47  shows  both  these  plots.  Compare  the  240  ksi  case  with  Figure  48  below. 
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Figure  47.  Membrane  Deflection  (Iiienl  =  239  and  240  ksi) 

The  change  in  E  seems  to  determine  when  one  layer  dominates  over  the  other  one. 

Getting  far  enough  from  the  “flip”  value,  it  is  noticed  the  E  value  changes  are  not 
affecting  the  final  solution  drastically.  This  can  be  seen  in  the  Figure  48,  which  has 
E=100  psi  on  the  left  and  E=10  ksi  on  the  right.  Observe  that  the  two  orders  of  magnitude 
difference  in  the  modulus  only  lead  to  a  deflection  change  of  a  small  fraction  of  a 
wavelength. 
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Figure  48.  Membrane  Deflection  (E^,em  =  100  psi  and  10  ksi) 

This  phenomenon  becomes  more  prominent  when  the  PVDF  modulus  was 
increased  to  2.8  Msi  and  2.8  Gsi.  Until  the  modulus  of  the  membrane  layer  was 
increased  to  similar  orders,  the  solutions  remained  nearly  the  same  for  all  membrane 
moduli,  as  seen  in  Figure  49,  plots  for  E  =  100,  10000  and  100000  from  left  to  right  for 
2.8  Msi. 
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Figure  49.  Membrane  Deflections  (Emem  =  100, 10000  and  100000  psi) 
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Once  the  membrane  modulus  was  increased  to  a  greater  order  of  magnitude  as  the 
PVDF  for  the  2.8  Msi  case,  the  solution  began  to  show  different  deflections,  but  they 
only  increased  slightly  (a  couple  of  wavelengths). 

One  last  part  of  this  test  is  to  show  the  effects  of  the  PVDF  layer  having  a  large 
value  of  E=2.8  Gsi.  The  three  plots  below  in  Figure  50  are  for  the  membrane  modulus  of 
1000,  106  and  109.  The  opposite  effect  as  that  of  the  other  PVDF  moduli  happened. 

When  the  membrane  modulus  was  increased  to  the  same  order  of  magnitude  as  the  PVDF 
layer,  the  deflection  actually  decreased.  The  scale  on  the  two  lower  order  moduli  is  104, 
while  it  is  10  for  the  PVDF  and  membrane  moduli  at  10  .  This  shows  how  much  the 
PVDF  layer  can  control  the  whole  system  without  the  influence  of  the  membrane  layer 
having  that  much  effect. 
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Figure  50.  Membrane  Deflections  (EWm  =  1000, 100000  and  109  psi) 
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The  last  test  performed  is  the  element  size.  This  is  another  variable  that  affects 
the  e  factor.  The  size  is  controlled  by  number  of  elements  across  the  membrane 
centerline,  i.e.  the  more  elements,  the  smaller  the  element  size  and  therefore,  the  larger 
the  e  .  The  program  is  setup  in  such  a  way  that  only  certain  element  numbers  can  be 
setup  with  certain  configurations.  For  the  configuration  that  has  been  used  in  all  the  tests 
up  to  this  point,  approximately  3  cm  center  region,  an  odd  amount  of  element  numbers  is 
required.  It  is  also  to  be  noted  that  while  the  element  number  is  chosen  to  have  a  certain 
value,  once  the  Zernike  coefficients  were  found,  the  solution  plots  created  are  always 
based  on  0.1”  increments  to  make  the  solution  more  continuous  as  is  the  case  in  the 
experimental.  This  should  be  an  easy  adjustment  if  it  is  desired  to  look  at  a  straight  line 
solution  between  nodes. 

This  test  performed  results  for  five  sets  of  element  numbers,  namely  3,  9,  15,  21, 
and  35  elements.  The  number  of  elements  across  the  whole  diameter  would  actually  be 
double  across  a  diameter,  but  since  this  was  an  axisymmetric  solution,  the  solution  covers 
across  one  radius.  While  many  other  tests  could  be  completed  in  this  arena,  these  five 
show  enough  of  the  solution  to  give  understanding  into  how  the  element  number  controls 
the  solution  by  completing  a  whole  cycle  of  deflection  plots,  as  will  be  described  below. 

The  first  test  had  the  lowest  amount  of  elements  possible  for  this  configuration, 
three.  During  the  process  of  understanding  and  learning  operation  of  the  program  it 
quickly  became  obvious  that  three  elements  were  going  to  be  a  problem.  A 
representative  test  using  the  following  properties,  I}mem=100  psi,  Efo/DF=261  ksi, 
tmem=0.006”,  tpvDF=0.003”,  N=1000,  300  V  in  center,  v  =0.3  for  both,  leads  produced  a 
solution  as  seen  in  Figure  5 1 . 
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Figure  51.  Membrane  Deflection  (3  Elements) 


Notice  the  very  large  deflection  with  this  amount  of  elements.  The  8  for  this  case  would 
be  very  small,  so  there  must  be  something  else  causing  the  error.  One  suggestion  for  this 
may  be  the  lack  of  rigidity  in  the  element  caused  by  being  too  large.  There  are 
configurations  that  the  three  elements  may  provide  reasonable  results,  but  at  this  time 
they  are  not  explored. 

For  the  time  being,  the  nine  element  case  is  going  to  be  skipped  for  sake  of  seeing 
what  happens  as  the  number  of  elements  is  increased.  The  15  element  case  is  therefore 
next.  The  same  material  properties  as  above  are  used  and  the  results  seem  reasonable  for 
the  low  modulus  of  the  membrane,  see  Figure  52.  It  shows  about  6  X  (approximately  4 
pm)  of  deflection  at  the  center.  This  is  more  than  is  desired,  but  the  solution  doesn’t 
seem  far  off.  This  is  the  beginning  of  the  divergence  on  the  high  number  of  elements. 
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Figure  52.  Membrane  Deflection  (15  Elements) 


Next,  21  elements  are  tested  and  the  deflection  is  five  times  greater  than  the  15 
element  solution,  see  Figure  53.  The  solution  illustrate  deflections  that  wouldn’t  be 
expected  for  the  relatively  low  voltage  conditions.  One  more  test  is  required  to  verify 
this. 
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Figure  53.  Membrane  Deflection  (21  Elements) 


The  number  of  elements  was  increased  by  a  little  more  than  50  percent  to  35 
elements.  Figure  54  illustrates  the  solution  ‘blowing  up’  again  for  the  high  number  of 
elements  having  nearly  105  X  of  displacement  or  on  the  order  of  14”,  which  is  clearly 
impossible  for  the  conditions  given  This  is  caused  by  the  elements  causing  a  large 
epsilon  Since  this  epsilon  term  is  fed  into  a  few  of  the  non-dimensional  terms,  which  are 
used  in  the  nonlinear  finite  element  code,  the  effects  of  it  increasing  in  the  integral  leads 
to  the  total  solution  becoming  invalid.  Further  study  of  this  effect  is  desired  and  required. 
The  reasoning  behind  the  element  rigidness  is  the  8  term.  As  the  element  size  becomes 
so  small,  e  becomes  very  large.  The  perturbation  solutionis  thereby  not  satisfied  by  the 
order  e  "  assumed.  This  is  a  local  effected  e  that  becomes  large  and  not  the  global  one; 
however,  it  must  still  be  satisfied  for  the  solution  to  remain  reasonable.  Epsilon  is 
directly  defined  by  the  thickness  and  length  of  each  element.  The  adjustment  of  this  term 
is  done  through  changing  the  values  of  these  two  parameters.  One  other  feature  in  this 
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study  is  the  “flip”  happening  from  negative  to  positive  to  negative  deflections.  The  two 
negative  solutions  are  discounted  as  poor  solutions  and  therefore,  not  considered  viable 
solutions. 
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Figure  54.  Membrane  Deflection  (35  Elements) 


Retreating  back  to  the  9-element  case,  it  is  seen  that  this  supplies  the  best  solution 
for  this  configuration.  The  deflection  plot  can  be  seen  in  Figure  55.  Notice  that  there  are 
just  a  few  wavelengths  of  deflection.  This  is  much  more  comparable  to  what  is  expected, 
on  the  order  of  a  micron  of  deflection. 
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Figure  55.  Membrane  Deflection  (9  Elements) 


The  above  study  dealing  with  the  element  size  shows  how  important  it  is  to  have 
elements  that  are  just  “rigid”  enough  to  provide  solutions  that  make  sense  and  are 
reasonable.  In  some  of  the  cases,  MATLAB  provided  warning  message  regarding  the 
size  vector  not  being  a  row  vector  with  integer  elements.  This  happened  at  35  elements 
and  could  have  caused  the  large  deflection  solution.  This  is  usually  caused  by  infinite 
values  in  the  vectors  and  is  just  another  verification  to  avoid  the  small  element  size. 

The  element  size  study  completes  the  parametric  study.  The  study  demonstrated 
how  there  is  probably  an  anomaly  in  the  code  when  it  comes  to  very  specific  values  of 
the  tension  and  modulus.  Around  these  values  the  deflection  reverses  signs.  Modulus 
and  tension  values  may  provide  poor  solutions  if  not  properly  applied.  Once  away  from 
the  anomalous  values,  the  solution  seemed  to  behave  and  actually  showed  what  would  be 
expected  as  the  tension  or  modulus  was  increased  or  decreased.  The  thickness  on  the 


other  hand  played  a  crucial  role  in  terms  of  the  magnitude  of  displacement,  granted,  it 
was  varied  on  the  same  scale  of  the  modulus  and  tension.  The  thicker  membrane  was  not 
a  valid  model,  which  also  happened  with  the  element  size  (or  element  number).  The 
increase  in  both  thickness  and  number  of  elements  caused  the  solution  to  grow  very 
rapidly.  There  are  obviously  values  for  all  the  variables  that  would  provide  excellent 
solutions  and  the  above  study  showed  the  sensitivity  to  changes  in  some  of  them. 

The  total  study  shows  the  influence  of  the  different  variables  in  the  final  solution. 
This  study  was  by  no  means  exhaustive,  but  it  covered  most  of  the  important  components 
required  before  moving  into  the  comparison  between  the  finite  solution  and  experimental 
data. 

Experimental  Test  Runs  Available 

Before  proceeding  directly  into  comparison  of  the  two  solutions,  a  short 

explanation  concerning  the  experimental  data  that  would  prove  useful  to  look  at  is 
included  here.  There  was  quite  a  bit  of  data  collected  by  Sobers  during  the  experimental 
stage.  As  described  in  Chapter  2,  there  were  two  “membrane”  test  articles  and  two  “stiff’ 
articles.  The  finite  element  code’s  capabilities  automatically  eliminated  the  two  stiff  test 
articles.  Ml  was  designed  mostly  to  learn  more  about  the  manufacturing  intricacies.  It 
could  be  used  as  a  test  article  to  compare  with  the  finite  element  code,  but  due  to  some  of 
the  concerns  Sobers  had  with  leak  over  of  the  piezo  layer  and  other  issues,  it  was 
determined  that  it  might  be  best  just  to  stick  with  M2.  M2  had  better  characteristics  in 
the  layout  of  the  piezo  zones  and  provided  more  data  in  the  axisymmetric  cases;  therefore, 
this  was  to  be  the  sole  test  article. 
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Sobers  had  performed  basically  four  tests,  not  including  the  OV  cases,  that  can  be 
easily  compared  (perfectly  axisymmetric).  Each  test  corresponded  to  a  different  voltage 
condition  in  300  V  increments  from  -600  V  to  600  V.  In  each  case,  after  the  positive  or 
negative  voltages  were  complete,  a  new  baseline  measurement  was  completed  at  0V. 

This  information  will  be  included  in  the  plots  below;  however,  the  difference  in  the  0V 
finite  element  solution  and  the  experimental  will  always  be  equal  to  the  experimental. 

This  is  due  to  the  fact  that  the  finite  element  solution  will  always  provide  no  deflection 
for  a  0  V  input  on  the  piezo  layer,  since  it  basically  represents  an  ideal  case.  The  plots 
provided  here  have  been  produced  using  Zernike  coefficients  and  polynomials  as 
explained  before.  The  boundary  is  actually  represented  at  x=0  and  the  center  at  x=3.  The 
y-axis  represents  deflection  or,  more  appropriate,  the  wavefront  deviation  from  flat.  This 
is  mentioned  since  there  are  times  when  the  experimental  boundary  is  not  located  at  the 
zero  mark  on  the  y-axis. 

The  order  that  appears  here  will  be  the  same  order  Sobers  used  in  the  testing 
procedure,  or  0V,  300V,  600V,  0V,  -300V,  600V,  0V.  The  contour  and  fully  3-D  plots 
from  the  experiment  for  both  the  full  Zernike  plots  and  the  reduced  data  will  be  included 
to  have  a  quick  reference.  The  finite  element  code  will  be  run  with  material  properties 
that  are  most  realistic  and  not  the  ones  that  produce  the  closest  solution.  The  properties 
used  for  this  solution  is  150  psi,  0.497,  and  0.25”  for  the  membrane  modulus,  Poisson’s 
ratio  and  thickness,  respectfully  and  261  ksi,  0.35  and  .002”for  the  same  PVDF  properties. 
Therefore,  there  will  be  plots  of  the  finite  element  code  that  potentially  could  be  orders  of 
magnitude  different  than  the  experimental.  If  this  is  the  case,  another  plot  of  the 
experimental  will  be  included,  so  the  shape  can  be  seen  more  clearly.  Other  solutions 
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may  be  included  representing  solutions  that  are  nearly  the  experimental.  In  these  cases, 
the  material  properties  will  be  given  that  provides  the  closest  solution. 

Experimental  and  Finite  Element  Comparison 

As  stated  before,  the  first  case  to  be  compared  (or  just  to  set  up  a  baseline)  is  the 
OV  case.  The  OV  experimental  plots  of  the  full  Zernike  and  reduced  Zemike,  along  with 
the  raw  data  are  included  in  Figure  56-Figure  57.  This  is  the  first  OV  case  before  any 
activation  was  applied  to  any  section  of  the  membrane. 


Figure  56.  Experimental  Contour  Plots  (Non -axisymme trie  and  Symmetric) 


Figure  57.  Experimental  Surface  Plots  (Non-axisymmetric  and  Symmetric) 
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The  above  plots  are  supplied  to  show  the  3-dimensional  and  contour  shapes,  as  well  as 
both  the  symmetric  and  non- symmetric  cases.  It  is  seen  that  there  is  some  difference  in 
wavefront  even  after  there  was  a  self- referencing  of  the  data.  This  could  be  shown  to  be 
“flat”  if  the  scale  on  the  z-axis  were  changed;  however,  it’s  important  to  show  that  there 
is  some  wavefront  error  on  some  level  even  if  it  is  very  small. 

The  following  plots  are  self-explanatory.  Figure  58  is  both  the  finite  element 
solution  defined  by  flat  line  at  zero  deflection  and  the  experimental  as  the  curved  line. 
Since  there  is  no  voltage  (OV)  input,  the  finite  element  solution  would  have  no  deflection 
as  represented  by  this  plot.  The  experimental  plot  is  actually  caused  by  the  difference  in 
the  wavefront  as  the  light  reflected,  and  therefore  can  have  a  non- flat  shape  for  no  voltage. 
The  plot  directly  below  in  the  same  figure  is  the  difference  between  the  two  plots. 
Obviously,  the  difference  in  this  case  is  equivalent  to  the  experimental  solution. 


Experimental  and  FEM 


Figure  58.  Reduced  Experimental  and  Finite  Solution  with  Difference 
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Figure  59  represents  the  total  solutions  of  the  35  Zemike  coefficients  for  the 
angles  defined  before  along  with  the  31  Zemike  coefficients,  these  two  plots  reduced  to  a 
single  line  and  the  modified  versions  where  the  deflections  corresponding  to  the 
centerlines  on  the  negative  x-axis  are  multiplied  by  -1,  defined  before.  These  plots  are 
based  solely  on  the  experimental  solution  and  are  included  as  a  reference  to  show  the 
effects  or  lack  there  of  the  non-axisymmetric  components.  The  scales  here  are  on  the 
order  of  0.2  X  . 
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Figure  59.  Non- Axisymme trie  Solutions 
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The  zero  volt  case  is  provided  not  to  try  to  verify  the  results  in  this  case,  but  so 
the  residual  effects  can  be  seen  after  the  voltage  has  been  applied.  It  also  presents 
information  about  the  initial  wavefront  before  testing.  The  next  two  tests  will  show  the 
program’s  capabilities  in  regards  to  a  voltage  input. 

The  first  voltage  condition  to  be  evaluated  is  a  positive  300  V.  For  now,  since  the 
material  properties  are  set  by  realistic  values  defined  above,  the  only  variable  that  could 
change  is  the  tension  term.  The  first  test  is  completed  with  a  tension  of  1000.  The  finite 
element  solution  produced  is  extremely  large  and  not  representative  of  the  actual  or 
experimental  solution.  Figure  60shows  the  two  solutions  together. 


x  10 


Experimental  and  FEM 
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Figure  60. 300V  Both  Solutions  and  difference 


Notice  that  the  finite  solution  dominates  the  plots.  Since  this  is  the  case,  the  plots  are 
separated  out  to  see  them  next  to  one  another,  as  seen  in  Figure  61. 
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Finrte  Element  Solution  Voltage  -  300 
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Figure  61.  FEM  and  Experimental  on  Separate  Plots 


The  next  step  is  to  see  if  there  is  a  possibility  that  a  reasonable  tension  can  be  found  that 
would  produce  a  solution  representative  of  the  experimental.  If  this  is  possible,  this 
tension  can  be  checked  against  all  the  other  voltage  conditions.  Unfortunately,  the 
tension  alone  could  not  bring  the  solution  to  the  same  order  without  cranking  it  up  to  an 

n 

unreasonable  value.  When  the  value  was  increased  to  8x10  ,  the  deflection  at  the  center 
was  almost  dead  on,  approximately  0.1  wavelengths  difference,  see  Figure  62. 
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Figure  62.  Both  Solutions  (T=8xl07) 


Following  the  order  Sobers  used  in  testing  the  membrane,  600  V  was  tested  next. 
The  large  tension  value  from  above  was  used  in  this  case.  The  solution  wasn’t  as  close, 
but  it  was  still  within  one  wavelength  of  difference  and  the  deflection  was  in  the  same 
direction.  This  can  be  seen  in  Figure  63. 
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Figure  63. 600V  Membrane  Deflection 


Tweaking  the  tension,  a  closer  approximation  can  be  found.  It  was  eventually 

o 

discovered  that  1x10  was  the  value  that  did  this.  However,  this  was  done  through  the 
use  of  knowing  the  experimental  solution  and  the  solution  was  still  about  1/2  X  difference. 
There  will  be  tests  that  there  is  no  experimental  solution  to  compare  and  the  tension  value 
cannot  be  adjusted  until  the  solution  matches;  therefore,  for  the  next  two  tests,  the  tension 
value  is  going  to  be  reduced  back  to  the  8x10  value. 

The  next  two  and  last  tests  for  the  comparison  with  M2  are  the  -300  and  -600  V 
cases.  The  tension  case  used  actually  produced  results  expressive  of  the  experimental. 
There  is  approximately  0.2  X  difference  between  the  two  solutions  for  -300V  and  nearly 
Vi  X  for  the  -600  V  case.  Both  of  these  plots  are  in  Figure  64. 
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Experimental  and  FEM 
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Figure  64.  Membrane  Deflection  (-300  and  -600V) 


The  above  tests  demonstrate  the  capabilities  of  the  finite  element  program  both  as 
a  separate  entity  and  when  in  comparison  with  an  actual  experimental  solution.  There  are 
many  other  tests  and  conditions  that  would  be  interesting  to  investigate.  Some  of  these 
will  be  mentioned  in  Chapter  5. 

On  the  whole,  the  finite  element  code  is  representative  of  the  experimental 
solution  for  a  reduced  Zemike  set.  It  was  possible  to  obtain  solutions  that  were  similar  to 
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the  experimental  solution;  however,  it  required  tweaking  of  a  variable,  in  this  case,  the 
boundary  tension  until  a  representative  value  was  obtained  to  match  the  first  voltage 
condition.  This  tension  value  produced  solutions  for  all  voltage  conditions  that 
represented  the  experimental  results.  As  the  parametric  study  showed,  the  many 
variables  available  can  have  drastic  effects  on  the  final  solution,  even  to  the  point  of 
reversing  the  direction  of  the  deflection.  For  this  reason,  it  is  important  to  have  the  best 
values  for  the  material  properties,  whether  this  means  actual  real-world  properties  or 
properties  of  a  nature  that  produces  close  solutions.  In  either  case,  if  there  is  not  an 
experiment  to  compare  with,  coming  back  to  these  tests  to  check  for  the  variable  values 
would  prove  useful  and  necessary,  so  as  to  avoid  creating  solutions  that  wouldn’t 
represent  a  real  experiment. 
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V:  Conclusions  and  Recommendations 


Finite  element  modeling  of  polymer  membrane-PVDF  structures  is  becoming 
more  necessary  with  the  technology  advancements  coming  down  the  pike.  This  study 
showed  the  beginnings  of  this  capability  in  symmetric  solutions  using  perturbation 
techniques.  In  many  cases,  a  symmetric  solution  may  be  enough  for  what  is  required  of 
the  mission,  implying  that  the  activation  zones  are  needed  to  change  the  global  shape  in  a 
symmetric  pattern.  If  this  is  the  case,  this  program  is  invaluable;  however,  there  may  be 
reason  to  have  non-axisymmetric  activation  zones  placed  in  different  regions  of  the 
mirror.  If  this  was  required  and  modeling  of  the  system  was  required  before  hand  to  see 
the  wavefront  shapes  possible,  the  symmetric  equations  would  not  provide  the  necessary 
conditions  and  would  therefore  need  to  be  modified.  This  may  be  accomplished  through 
extending  formulation  to  two  dimensions  and  consequently  the  shape  functions  would  be 
better  suited  to  the  solution.  If  this  new  program  was  developed,  it  would  be  interesting 
to  recheck  the  experimental  solutions  without  reducing  them  to  the  symmetric  case. 

When  the  program  was  first  developed  and  the  experimental  tests  were  being 
completed,  it  was  known  that  they  were  related  projects  and  it  was  determined  to  have 
axisymmetric  experiments  completed.  Now  that  more  work  has  been  completed  with  the 
program,  more  experimental  tests  with  different  PVDF  configurations  and  voltage 
conditions  would  supply  more  information  regarding  the  capability  of  the  program  and 
give  more  to  compare  against.  It  would  also  show  experimentally  what  may  be  possible 
in  shape  control  on  the  optical  level. 
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If  this  were  completed,  different  configurations  in  the  finite  element  solution 
could  be  investigated,  rather  than  just  the  center  circle.  There  are  many  different 
symmetric  configurations  of  concentric  circles  around  the  membrane  that  could  produce 
solutions.  If  this  study  is  completed,  it  would  be  beneficial  to  rewrite  part  of  the  program 
to  have  the  ability  to  have  different  voltage  conditions  in  each  of  the  activation  zones.  As 
it  stands  now,  each  zone  that  is  given  activation  (see  Appendix  B)  receives  the  voltage 
defined  at  the  beginning.  In  a  real  life  situation,  there  may  be  a  few  voltage  sources  that 
supply  different  activation  zones  and  therefore  would  prove  useful  to  be  able  to  have  this 
in  the  program  as  well. 

The  parametric  study  usually  looked  at  changing  one  variable  at  a  time  to 
investigate  the  singular  effects  of  a  variable.  It  may  be  interesting  to  investigate  various 
changes  together;  so  as  to  optimize  the  material  properties  to  match  the  experimental 
results.  The  results  of  the  comparison  showed  decent  results,  but  it  was  not  looked  at 
how  to  match  the  solution  exactly.  It  may  be  beneficial  to  write  a  short  script  to  converge 
on  the  closest  zero  difference  solution  possible  by  varying  the  modulus,  tension,  or  any  of 
the  other  material  properties  together.  The  solution  to  this  study  would  be  the  material 
and  system  properties  that  gave  the  best  solution  for  each  voltage  condition.  They  could 
be  compared  to  see  if  the  properties  were  in  the  same  ballpark  for  each  voltage.  This 
should  be  accomplished  by  optimizing  over  all  voltage  experimental  results 
simultaneously. 

As  shown,  the  solution  can  be  very  sensitive  to  increasing  £  .  For  this  reason,  a 
large  tension  was  required  for  the  actual  tests  when  the  membrane  was  *4”  thick.  The 
viscosity  and  curing  process  of  the  RTV  615  is  the  reason  for  the  membrane  being  that 
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thick.  Usually,  the  membranes  are  more  on  the  order  of  tens  of  microns.  These  polymer 
membranes  (one  such  is  SRS  Technologies,  Inc.  CP-1)  available  are  usually  of  better 
surface  quality  than  the  RTV  615  and  therefore  some  of  the  non-axisymmetric  wavefront 
deviations  may  be  reduced  or  eliminated  if  they  were  used.  CP- 1  or  other  material  may 
prove  useful  to  test  in  the  same  manner  of  RTV  615.  If  nothing  else,  it  would  provide 
more  data  that  could  be  compared.  It  could  also  see  if  CP- 1  is  a  more  viable  option  for 
wavefront  shape  control  or  if  it  would  prove  to  be  as  responsive  to  the  PVDF  changes  as 
the  RTV615.  The  material  used  here  is  just  as  important  as  the  ability  to  control  and 
model  it. 

The  study  completed  above  demonstrated  the  basic  capabilities  of  the  finite 
element  code.  There  are  obviously  many  shortcomings  in  the  solutions.  The  first  test 
completed  was  to  discover  the  tension  effects  on  the  membrane.  The  only  solutions 
similar  to  the  experimental  results  were  when  the  tension  was  increased  to  infeasible 
values.  When  this  value  was  used  under  other  voltage  conditions,  the  solutions  seemed  to 
remain  representative  of  the  experimentally.  Although  the  number  used  to  obtain  this 
solution  is  unreasonable,  it  may  be  useful  to  continue  to  use  this  value  whenever  an  RTV 
membrane  is  tested.  A  different  membrane  material  may  require  a  different  tension  value 
and  even  a  different  RTV  membrane  may  have  a  different  value.  A  method  to  find  a 
value  experimental  could  prove  useful.  A  vibration  test  may  be  one  way  this  could  be 
accomplished.  For  comparison  to  the  experimental,  this  was  the  only  variable  varied, 
since  the  material  properties  were  set  by  the  values  of  Chapter  2  Material  Testing. 

The  parametric  study  showed  how  changes  in  different  properties  sometimes  lead 
to  extremely  large  changes  and  other  times  hardly  affected  the  solution.  The  changes  in 
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the  modulus  fell  into  the  latter  category.  There  was  an  area  of  moduli  that  had  large 
changes;  however,  once  away  from  these  values,  the  deflection  plots  are  nearly  the  same 
order.  As  with  many  of  the  properties,  there  was  a  value  that  had  a  flip  from  positive  to 
negative  displacement.  The  rationale  behind  this  is  unknown  and  must  correspond  to 
some  function  in  the  code  or  the  dominating  effect  of  one  layer  of  the  piezo -polymer 
layup. 

As  stated  before,  most  of  the  material  properties  have  this  displacement  flip. 
Although  written  above  for  the  specific  comparison,  the  tension  was  checked  out  in  the 
parametric  study  first.  As  it  was  increased  to  a  very  large  value,  it  was  shown  to  produce 
very  small  displacement  plots  (~10‘6  A, ).  It  also  had  a  flip  at  a  relatively  small  value, 
namely  3242.5.  On  either  side  of  this  value,  the  displacement  increased.  Values  less 
than  the  3200  value  the  displacement  kept  increasing  as  N  approached  zero.  The  other 
side  increased  for  a  while,  then  started  to  decrease  in  displacement  after  that  point,  until 
the  displacement  was  nearly  zero.  Away  from  the  anomaly  tension  value,  the  membrane 
behaves  as  it  should,  i.e.  low  tension  will  allow  greater  displacement  and  higher  boundary 
tension  will  produce  very  small  displacement.  The  anomaly  is  changed  from  the  3242.5 
value  if  the  modulus  is  increased  or  the  thickness  decreased. 

The  thickness  of  the  membrane  also  demonstrated  quite  a  bit  of  change  for  very 
small  changes.  The  one  case  demonstrated  that  increasing  the  thickness  from  0.0066”  to 
0.009”  lead  to  a  5000  times  increase  in  displacement.  This  sensitivity  in  thickness  is 
most  likely  explainable  in  terms  of  e  .  As  thickness  is  increased,  the  £  term  is  also 
increasing.  Once  it  becomes  larger  than  a  certain  size,  it  will  cause  the  solution  to  blow 
up.  This  is  supposedly  what  happened  in  this  case.  This  sensitivity  to  £  needs  to  be 
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further  investigated  to  possibly  find  a  way  to  reduce  its  effects.  The  £  sensitivity  was 
also  noticed  in  the  other  component  of  element  size,  the  number  of  elements  used. 

Up  to  this  point,  9  elements  were  used  for  all  the  tests.  However,  it  is  important 
to  see  what  happens  using  more  and  less  elements.  The  sensitivity  to  £  was  quickly  seen 
when  the  number  of  elements  was  increased;  thereby,  the  element  length  to  thickness  was 
decreased.  As  more  elements  were  investigated,  the  displacement  kept  increasing  to 
areas  that  were  completely  impossible  (105  X  ).  Nine  elements  always  produced  what 
seems  like  a  reasonable  solution.  Decreasing  the  number  of  elements  also  caused  a 
problem  in  the  code.  When  three  elements  were  used,  the  displacement  was  1014  X  . 

There  is  some  speculation  why  this  happened.  The  PVDF  layer  may  be  dominating,  but 
it  is  more  likely  dealing  with  the  element  size  being  so  large  that  the  effects  of  the 
activation  zone  is  not  able  to  produce  a  solution.  Another  possibility  may  reside  in  the 
non-dimensionalizing  of  the  terms  and  how  change  one  dimensioned  term  changes  a  few 
non- dimensioned  terms.  These  terms  are  solved  through  a  nonlinear  integral  and  could 
cause  many  of  the  problems  with  the  code  not  functioning  properely  for  a  few  of  the 
cases.  There  are  many  other  tests  in  this  study  that  could  be  complete;  however,  the 
program’s  capabilities  were  shown  in  it  and  the  effects  of  different  properties  investigated. 

Overall  this  study  was  enlightening  to  see  the  modeling  abilities  in  comparison  to 
experiments  and  as  a  stand  alone  program  in  the  parametric  study.  Since  this  is  still  a 
relatively  new  field  of  study  there  is  much  to  be  completed  in  the  future  and  the  above 
discussion  provides  a  starting  point,  but  doesn’t  cover  all  that  will  be  desired  in  further 
studies. 
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Appendix  A:  Finite  Element  Code 


The  following  appendix  is  the  different  parts  of  the  code  along  with  comments 
about  their  operation  or  purpose.  All  of  this  code  along  with  other  code  is  supplied  on  the 
program  CD. 

This  is  a  sample  of  the  experimental  data  that  produces  the  Zemike  plots.  There 
is  a  complimentary  raw  data  file;  however,  it  is  just  a  long  column  of  numbers,  so  it  is  not 
supplied  here.  Any  combination  of  Zemike  coefficients  could  be  typed  into  a  file  similar 
to  this  and  when  this  file  is  mn  in  MATLAB®,  will  create  4  plots  of  raw  data,  smoothed 
raw  data,  and  two  Zemike  plots. 

%  M2_Tlc_300v 
%  Mon  Dec  10  21:51:53  2001 
%  Zygo  Zemike  Coefficients 
%  Obscuration  Ratio  =  0.0000 
%  Index  Coefs(microns)  Equation 
%  function  [z]  =  threeOOcenter 
Z=[ 

1  0.000000  %rcos(t)(X  Tilt) 

2  0.000000  %  rsin(t)  (Y  Tilt) 

3  0.000000  %2rA2-l  (Focus) 

4  0.000000  %  rA2cos(2t)  (0  Astigmatism) 

5  0.000000  %  rA2sin(2t)  (45  Astigmatism) 

6  0.000000  %  (3rA2 -2 )rcos(t)  (X  Coma) 

7  0.000000  %  (3rA2-2)rsin(t)  (Y Coma) 

8  0.152906  %  6rA4-6rA2+l  (Spherical) 

9  0.000000  %  rA3cos(3t) 

10  0.000000  %  rA3sin(3t) 

11  0.000000  %(4rA2-3)rA2cos(2t) 

12  0.000000  %  (4rA2 -3 )rA2sin( 2t) 

13  0.000000  %  (10rA4-12rA2+3)rcos(t) 

14  0.000000  %  ( I0rA4-12rA2+3)rsin(t) 

15  0.188982  %  20rA6  -30rA4+ 12rA2  -1 

16  0.000000  %  rA4cos(4t) 

17  0.000000  %rA4sin(4t) 

18  0.000000  %  (5rA2-4)rA3cos(3t) 

19  0.000000  %  (5rA2-4)rA3sin(3t) 

20  0.000000  %  (15rA4-20rA2+6)rA2cos(2t) 

21  0.000000  %  ( 1 5 rA4 -20rA2+6)rA2sin(2t) 

22  0.000000  %  (35^6-60^4+30^2-4)1-^(1) 

23  0.000000  %(35rA6-60rA4+30rA2-4)rsin(t) 

24  0.142927  %  70rA8-140rA6+90rA4-20rA2+l 

25  0.000000  %  rA5cos(5t) 

26  0.000000  %  rA5sin(5t) 

27  0.000000  %(6rA2-5)rA4cos(4t ) 

28  0.000000  %  ( 6rA2 -5 )rA4sin(4t) 

29  0.000000  %  (21rA4-30rA2+10)rA3cos(3t) 

30  0.000000  %(21rA4-30rA2+10)rA3sin(3t) 

31  0.000000  %  (56rA6-105rA4+60rA2-10)rA2cos(2t) 

32  0.000000  %  (56rA6-105rA4+60rA2-10)rA2sin(2t) 

33  0.000000  %  ( 126rA8 -280rA6+ 2 1 0rA4 -60rA2 +5)rcos(t) 
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34  0.000000  %  ( 126rA8 -280^6+210^4 -60r*2+5)rsin( t) 

35  0.123484  %  252^10-630^8+560^6-210^4+30^2-1 

]: 


data = textread(  'M2_T 1  c_3 OOv. dat',  '%f ) ; 

%factor(length(Tlg))  %  [3331  37] 

rows- 111; 

cols=93; 

zsize=94; 

smoothlevel=l ; 

pupil=0; 

shift=[0  0]: 

MLM=0.133; 

ratio=7.6327; 

trim=l; 

makeplots  %  Command  that  calls  script  file  to  plot  surface  plots 
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This  file  takes  the  information  from  the  raw  data  and  Zernike  file  above  and  creates  the 
plots.  It  is  slightly  modified  from  its  original  form  to  put  all  the  plots  on  the  same  figure. 
If  surface  plots  are  desired,  contour  should  be  replaced  with  surf. 

%  makeplots.m  subroutine 

%  Copyright  2002  -  Michael  Sobers  -  All  rights  reserved. 

%  Modified  by  Brian  Lutz  2004 
z=z(:,2); 

for  j  =  l:cols 

data_surf(  1  :rows,j  )=data(  (j  -1  )*rows+ 1  :j  *rows); 
end 


if  trim==l 

temp=zeros(rows+2,cols+2); 

temp(2:rows+l,2:cols+l)=data_surf; 

data_surf=temp; 

while  (sum(data_surf(2,:))==0&sum(diff(data_surf(2,:)))==0) 
data_surf(2:size(data_surf,l)-l,:)=data_surf(3:size(data_surf,l),:); 
end 
r=2; 

while  (sum(data_surf(r,:))~=Olsum(diff(data_surf(r,:)))~=0) 
r=r+l; 
end 

data_surf=data_surf(  1  :r,:); 

while  (sum(data_surf(:,2))== 0&sum(diff(data_surf(:,2)))==0) 
data_surf(  :,2:  size(data_surf,2)  - 1  )=data_surf(:  ,3 :  size(data_surf,2)); 
end 
c=2; 

while  (sum(data_surf(:,c))~=Olsum(diff(data_surf(:,c)))~=0) 
c=c+l; 
end 

data_surf=data_surf  ( : ,  1 :  c) ; 
rows=size(data_surf,l); 
cols=size(data_surf,2); 
end 

data_surf=data_surf./(0.633); 

shift=round(shift/(MLM*ratio)); 

zsurface  =  myzern(z,[rows  cols],  pupil,  shift); 

array=zsurface; 

array=array./(0.633); 


if  mkpts  ==  1 
figure 

subplot(2,2,l),  contour] [0:size(array,2)-l]*MLM*ratio,[0:size(array,l)-l]*MLM*ratio,array); 

view([-19,68]) 

axis([0  110  0  110  -1  1.5]) 

colormap(copper) 

colorbar 

shading  interp 
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xlabel=(’Surface  Location  -  mm'); 
ylabel=(’Surface  Location  -  mm'); 
zlabel=(’Surface  Height  -  Wavelengths  (633nm)'); 
end 

%for  j  =  lxols 

%  data_surf_flip(:,j)=data_surf(:,cols+l  -j); 

%end 

%for  r  =  l:rows 

%  data_surf_rotated(:,r)  =  data_surf_flip(rows+l-r,:)'; 
%end 


if  mkpts  ==  1 

subplot(2,2,2),  contour( [Oxols  -1]  *MLM*ratio, [0:rows - 1]  *MLM*ratio,data_surf); 
view([-19,68]) 

%axis([0  110  0  110-1  1.5]) 
colormap(copper) 
colorbar 
shading  interp 
end 

smooth_data  =  interpolate2(data_surf,smoothlevel); 
if  mkpts  ==  1 

subplot(2, 2, 3), contour!  [Oxols -l]*MLM*ratio,[0:rows-l]*MLM*ratio,smooth_data); 
view([-19,68]) 

%axis([0  110  0  110-1  1.5]) 
colormap(copper) 
colorbar 
shading  interp 
end 

mask=zeros(rows,  cols); 

for  r=l:  rows 
for  c=  lxols 
if  smooth_data(r,c)~=0 
mask(r,c)=l; 
end 
end 
end 

%  figure(3) 

H=axis; 


if  mkpts  ==  1 

subplot(2, 2, 4), contour!  [0:size(array,2>l]*MLM*ratio,[0:size(array,l)-l]*MLM*ratio, array.  *mask); 

view([-19,68]) 

axis(H); 

%axis([0  110  0  110-2  2]) 
colormap(copper) 
colorbar 
shading  interp 
end 
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This  is  the  first  of  the  finite  element  code  files.  It  contains  the  material  properties  and  it 
prepares  them  by  placing  them  in  matrix  form  before  the  information  is  exported  to 
another  file  that  runs  the  calculations. 

function  [v,x,d,eps]=testr2(numel,p,t,v, order, verb, dyn.mode, alpha, beta, config); 

innerR=0.0001; 

outerR=3.0001; 

%  youngsperlayer=diag([406000  261000  261000]); 

%  thickperlayer=diag( [0.006  0.003  0.003]); 

%  poissonperlayer=diag([0.3  0.3  0.3]); 

%  densityperlayer=diag([0.1  0.12  0.12]); 

%  d31perlayer=diag([0  4  4]*10A(-8)); 

%  thick=thickperlayer*ones(3,numel); 

%  dens=densityperlayer*ones(3,numel); 

%  d31=d31perlayer*ones(3,numel); 

%  youngs=youngsperlayer*ones(3,numel); 

%  poisson=poissonperlayer*ones(3,numel); 
youngsperlayer=diag(  [  1 00  261000]); 
thickperlayer=diag([.006  .003]); 
poissonperlayer=diag([.3  .3]); 
densityperlayer=diag([.l  .12]); 
d31perlayer=diag([0  4]*10A(-8)); 
nodes=innerR+[0: 1  :numel]  *(outerR-innerR)/numel; 
thick=thickperlayer*ones(2,numel); 
dens=densityperlayer*ones(2,numel); 
d3  l=d3  lperlayer*ones(2,numel); 
youngs=youngsperlayer*ones(2,numel); 
poisson=poissonperlayer*ones(2,numel); 
tens=t/sum(sum(thickperlayer)); 

if  mod(numel,2)  ==  0 
if  config  ==  1 

volts=v*[ones(2,numel-l)  zeros(2,l)];  %Configuration  1 
elseif  config  ==  0 

volts=v*[ones(2,numel)];  %Configuration  0 
end 

elseif  mod(numel,2)  >  0 
if  config  ==  4 

volts=v*[zeros(2,numel/3)  ones(2,numel/3)  zeros(2,numel/3)];  %  Configuration  4 
elseif  config  ==  3 

volts=v*[-ones(2,numel/3)  ones(2,numel/3)  zeros(2,numel/3)];  %Configuration  3 
elseif  config  ==  2 

volts=v*[ones(2,numel/3)  ones(2,numel/3)  zeros(2,numel/3)];  %Configuration  2 
elseif  config  ==  5 

volts=v*[ones(2,numel/3)  zeros(2,numel/3)  zeros(2,numel/3)  zeros(2,numel/3)  zeros(2,numel/3)];  % 
3.048  cm  circle  in  center 
end 
end 

[v,x,d,eps]=memfemr2(order,p,tens,nodes,thick,dens,youngs,poisson,  d31,  volts, verb.dyn.mode, 3, alpha, 
beta); 
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The  following  file  is  where  the  calculations  are  performed.  No  modifications  involved 
with  the  simple  problem  should  be  made  to  this  file.  There  are  a  few  files  run  by  this  file 
and  are  bcated  in  the  Private  folder.  This  file  and  the  corresponding  files  contained 
within  are  where  the  theoretical  development  was  programmed. 

function  [vr.xr,  dr,eps]=memfemr2(order,  p,  tens,  nodes,  thick,  dens,  youngs,  poisson,  d31,  volts,  verbose, 
dyn,  mode,  mx,  alpha,  beta) 

%  bs3fem  -  Finite  Element  simulation  of  plate- membrane 
%  written  by  :  Capt  James  Rogers,  AFIT,  USAF 

% 

%  order  =  order  of  perturbation  levels  terms  to  include  (0  to  3) 

%  p  =  pressure  differential  per  unit  length  applied  to  beam 
%  tens  =  pretension  applied  to  beam 
%  nodes  =  distance  from  left  end  of  each  element  interface 
%  thick  =  thickness  of  layers 
%  dens  =  layer  material  density  (mass/LA3) 

%  youngs  =  layer  Young's  Modulus  (F/LA2) 

%  poisson  =  layer  Poisson's  Ratio 
%  d31  =  layer  piezo  coefficients 

%  volts  =  layer  voltage  potential 
%  verbose  =  0:  No  Plots  Produced 
%  1:  All  Plots  Created 

%  dyn  =  0:  Static  Only 
%  1 :  Dynamic  V  V  V  V 

%  mode  =  0  :  Clamped  Natural  Modes  Returned 
%  >0:  assumed  center  mode 

%  mx  =  maximum  modes  returned 
% 

%  vr  =  Eigenvectors  (Dynamic) 

%  Analytic/Total  Static  FEM/Only  1st  Level  FEM  Solutions  (Static) 

%  dr  =  Eigenvalues  (Dynamic) 

%  =  Corresponding  Axial  Locations  (Static) 

% 


num=length(nodes)-l; 
len  =  nodes(length(nodes)); 
if  max(size(thick))  ==  num 
numlayers  =  min(size(thick)); 
else 

numlayers  =  max(size(thick)); 
end 

for  i=l:num 
hi=0; 
hbu=0; 
hbl=0; 

for  j=l:numlayers 

hbu=hbu+thick(j,i)*(2*hi+thick(j,i)); 

hbl=hbl+2*thick(j,i); 

hi=hi+thick(j,i); 

end 

hbar=hbu/hbl; 

%Dimensional  Values 
els(i).p  =p; 
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els(i).rl  =  nodes(i); 

els(i).r2  =  nodes(i+l); 

els(i).r  =  (nodes(i+l)+nodes(i))/2; 

els(i).len  =  (nodes(i+l)-nodes(i)); 

els(i).rhoa  =  0; 

els(i).t  =  0; 

els(i).ehll  =0; 

els(i).ehl2  =  0; 

els(i).ehe  =  0; 

els(i).ehzll  =0; 

els(i).ehzl2  =  0; 

els(i).ehze  =  0; 

els(i).dll  =0; 

els(i).dl2  =0; 

els(i).nz  =  0; 

els(i).n  =  0; 

hi=0; 

for  j = 1 :  numlayers 

zl  =  thick(j,i)*(2*hi+thick(j,i)-2*hbar)/2; 

z2  =  (thick(j,i)*(thick(j,i)A2+3*hi*thick(j,i)+3*hi*2)- 

3*hbar*thick(j,i)*(2*hi+thick(j,i))+3*hbarA2*thick(j,i))/3; 

els(i).rhoa  =  els(i).rhoa+dens(j,i)*thick(j,i); 
els(i).ehll  =  els(i).ehll  +  youngs(j,i)*thick(j,i)/(l-poisson(j,i)A2); 
els(i).ehl2  =  els(i).ehl2  +  youngs(j,i)*poisson(j,i)*thick(j,i)/(l-poisson(j,i)A2); 
els(i).ehe  =  els(i).ehe  +  youngs(j,i)*d31(j,i)*volts(j,i)/(l-poisson(j,i)); 
els(i).ehzll  =  els(i).ehzll  +  youngs(j,i)*zl/(l-poisson(j,i)A2); 
els(i).ehzl2  =  els(i).ehzl2  +  youngs(j,i)*poisson(j,i)*zl/(l-poisson(j,i)A2); 

els(i).ehze  =  els(i).ehze  +  youngs(j,i)*zl*d31(j,i)*volts(j,i)/thick(j,i)/(l-poisson(j,i)); 
els(i).dll  =els(i).dll  +  youngs(j,i)*z2/(l-poisson(j,i)A2); 

els(i).dl2  =  els(i).dl2  +  youngs(j,i)*poisson(j,i)*z2/(l-poisson(j,i)A2); 

els(i).nz  =  els(i).nz  +  tens*zl; 

els(i).n  =  els(i).n  +  tens*thick(j,i); 

els(i).t  =  els(i).t+thick(j,i); 
hi  =  hi+thick(j,i); 

end 

%Scaled  Values 

els(i).eta2  =  els(i).ehl  l/els(i).n; 

els(i).eta  =  sqrt(els(i).eta2); 

els(i).eps  =  sqrt(els(i).dll/els(i).ehll)/els(i).len; 

els(i).dl2s  =  els(i).dl2/els(i).ehl  l/els(i).epsA2*numA2; 
els(i).ps  =  p*els(i).len/els(i).ehl  l/els(i).epsA2*numA2; 

els(i).ehes  =  els(i).ehe*els(i).lenA2/els(i).ehl  l/els(i).epsA2; 
els(i).ehzl  Is  =  els(i).ehzll/els(i).ehl  l/els(i).len; 
els(i).ehzl2s  =  els(i).ehzl2/els(i).ehl  l/els(i).len; 
els(i).ehzes  =  els(i).ehze*els(i).lenA2/els(i).ehl  l/els(i).epsA3; 
els(i).nzs  =  els(i).nz/els(i).ehl  l/els(i).len; 
els(i).nubar  =  els(i).ehl2/els(i).ehl  1; 
eps(i)  =  els(i).eps; 

end 

disp=3;  %nodes  per  element 
dpn=2;  %displacements  per  node 

kgO=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 
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mgO=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 

pgO=zeros(dpn*((disp-l)*num+l),l); 

w01=zeros(dpn*((disp-l)*num+l),l); 

if  order  >  0 

kgl=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 
mgl=zeros(dpn*((disp-l)*nurrH-l),dpn*((disp-l)*num+l)); 
pgl=zeros(dpn*((disp-l)*num+l),l); 
wl  l=zeros(dpn*((disp-l)*num+l),l); 
if  order  >  1 

kg2=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 

mg2=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 

pg2=zeros(dpn*((disp-l)*num+l),l); 

w21=zeros(dpn*((disp-l)*num+l),l); 

if  order  >  2 

kg3=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 

mg3=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 

pg3=zeros(dpn*((disp-l)*num+l),l); 

w31=zeros(dpn*((disp-l)*num+l),l); 

end 

end 

end 

chunk=l; 

ds=zeros(chunk*num+ 1  ,order+ 1 ); 
x=zeros(chunk*num+l,l); 
ws=zeros  (chunk*num+ 1  ,order+ 1 ) ; 
w=zeros(chunk*num+l,l); 

for  i=l:num 

nodeids=[(disp-l)*i-(disp-2)  (disp-l)*i+(disp-2)  (disp-l)*i]; 
kO  =  KOfr(els(i).eps,els(i).r,els(i).len)  ... 

+  K0gr(els(i).eps,els(i).r,els(i).len)/els(i).eta2; 
rnO  =  MOr(els(i).eps,els(i).r,els(i).len); 
pO  =  els(i).ps*POr(els(i).eps,els(i).r,els(i).len); 
kgO  =  BuildStiffness(kgO,kO/els(i).len,nodeids); 
mgO  =  BuildMass(mgO,mO/els(i).len,nodeids); 
pgO  =  BuildForce(pgO,els(i).epsA2*pO,nodeids); 
if  order  >  0 

ml  =  Mlr(els(i).eps,els(i).r,els(i).len); 
kl  =  Klfr(els(i).eps,els(i).r,els(i).len)  ... 

+  Klgr(els(i).eps,els(i).r,els(i).len)/els(i).eta2... 

+  els(i).dl2s*Kldr(els(i).eps,els(i).r,els(i).len); 
fl  =  els(i).ehzes*FOr(els(i).eps,els(i).r,els(i).len); 
pi  =  els(i).ps*Plr(els(i).eps,els(i).r,els(i).len); 

kgl  =  BuildStiffness(kgl,els(i).eps*kl/els(i).len,nodeids); 
mgl  =  BuildMass(mgl,els(i).eps*ml/els(i).len,nodeids); 
pgl  =  BuildForce(pgl,els(i).epsA3*(fl-pl),nodeids); 
if  order  >  1 

m2  =  M2r(els(i).eps,els(i).r,els(i).len); 
k2  =  K2fr(els(i).eps,els(i).r,els(i).len) ... 

+  K2gr(els(i).eps,els(i).r,els(i).len)/els(i).eta2... 

+  K2nr(els(i).eps,els(i).r,els(i).len)... 

+  els(i).dl2s*K2dr(els(i).eps,els(i).r,els(i).len); 
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k2e  =  els(i).ehes*KOgr(els(i).eps,els(i).r,els(i).len); 

kg2  =  BuildStiffness(kg2,els(i).epsA2*(k2-k2e)/els(i).len,nodeids); 
mg2  =  BuildMass(mg2,els(i).epsA2*m2/els(i).len,nodeids); 
f2  =  els(i).ehzes*Flr(els(i).eps,els(i).r,els(i).len); 
p2  =  els(i).ps*P2r(els(i).eps,els(i).r,els(i).len); 

pg2  =  BuildForce(pg2,els(i).epsA4*(f2-p2),nodeids); 
if  order  >  2 

m3  =  M3r(els(i).eps,els(i).r,els(i).len); 
k3  =  K3fr(els(i).eps,els(i).r,els(i).len)  ... 

+  K3gr(els(i).eps,els(i).r,els(i).len)/els(i).eta2... 

+  K3nr(els(i).eps,els(i).r,els(i).len)... 

+  els(i).dl2s*K3dr(els(i).eps,els(i).r,els(i).len); 
k3e  =  els(i).ehes*Klgr(els(i).eps,els(i).r,els(i).len); 
kg3  =  BuildStiffness(kg3,els(i).epsA3*(k3  -k3e)/els(i).len,nodeids); 
mg3  =  BuildMass(mg3,els(i).epsA3*m2/els(i).len,nodeids); 

f3  =  els(i).ehzes*F2r(els(i).eps,els(i).r,els(i).len); 
p3  =  els(i).ps*P3r(els(i).eps,els(i).r,els(i).len); 

pg3  =  BuildForce(pg3,els(i).epsA5*(f3-p3),nodeids); 

end 

end 

end 

end 

%  clamped 

%  sort  Oth  order  matrices  for  undetermined  displacements 

mg01=mg0; 

templ=mg01(l,:); 

mg01(l,:)=mg01(2,:); 

mg01(2,:)=templ; 

temp2=mg01(:,l); 

mg01(:,l)=mg01(:,2); 

mg01(:,2)=temp2; 

kg01=kg0; 

templ=kg01(l,:); 

kg01(l,:)=kg01(2,:); 

kg01(2,:)=templ; 

temp2=kg01(:,l); 

kg01(:,l)=kg01(:,2); 

kg01(:,2)=temp2; 

Pg0(4)=p  g0(4)+pg0(2); 

pg0(  1  )=pg0(  l)+pg0(2)*els(  l).len/2; 

Pg01=pg0; 

temp=pg01(l); 

PgO  1  ( 1  )=pg0 1(2); 
pg01(2)=temp; 

m0c=mg01(dpn:dpn*(disp-l)*num,dpn:dpn*(disp-l)*num); 

k0c=kg01(dpn:dpn*(disp-l)*num,dpn:dpn*(disp-l)*num); 

p0c=pg01(dpn:dpn*(disp-l)*num); 

if  order  >  0 
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%  sort  1st  order  matrices  for  undetermined  displacements 

mgll=mgl; 

templ=mgl  1(1,:); 

mgll(l,:)=mgl  1(2,:); 

mgll(2,:)=templ; 

temp2=mgl  1(:,1); 

mgll(:,l)=mgll(:,2); 

mgll(:,2)=temp2; 

kgll=kgl; 
templ=kgl  1(1,:); 
kgl  1(1,  :)=kgl  1(2,:); 
kgl  l(2,:)=templ; 
temp2=kgl  1(:,1); 
kgll(:,l)=kgll(:,2); 
kgl  l(:,2)=temp2; 

pgll=pgl; 
temp=pgll(l); 
pgl  l(l)=pgl  1(2); 

Pgl  l(2)=temp; 

mlc=mgll(dpn:dpn*(disp-l)*num,dpn:dpn*(disp-l)*num); 
klc=kgll(dpn:dpn*(disp-l)*num,dpn:dpn*(disp-l)*num); 
plc=pgl  l(dpn:dpn*(disp-l)*num); 
clc  =  alpha*kOc  +  beta*mOc; 
if  order  >  1 

%  sort  1st  order  matrices  for  undetermined  displacements 

mg21=mg2; 

templ=mg21(l,:); 

mg21(l,:)=mg21(2,:); 

mg21(2,:)=templ; 

temp2=mg21(:,l); 

mg2 1  (: ,  1  )=mg2 1  ( :  ,2); 

mg21(:,2)=temp2; 

kg21=kg2; 

templ=kg21(l,:); 

kg21(l,:)=kg21(2,:); 

kg21(2,:)=templ; 

temp2=kg21(:,l); 

kg21(:,l)=kg21(:,2); 

kg21(:,2)=temp2; 

pg21=pgl; 

temp=pg21(l); 

pg21(l)=pg21(2); 

pg21(2)=temp; 

m2c=mg21(dpn:dpn*(disp-l)*num,dpn:dpn*(disp-l)*num); 

k2c=kg21(dpn:dpn*(disp-l)*num,dpn:dpn*(disp-l)*num); 

p2c=pg21(dpn:dpn*(disp-l)*num); 

c2c  =  alpha*klc  +  beta*mlc; 

if  order  >  2 
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%  sort  1st  order  matrices  for  undetermined  displacements 

mg31=mg3; 

templ=mg31(l,:); 

mg31(l,:)=mg31(2,:); 

mg31(2,:)=templ; 

temp2=mg21(:,l); 

mg31(:,l)=mg31(:,2); 

mg31(:,2)=temp2; 

kg31=kg3; 

templ=kg31(l,:); 

kg31(l,:)=kg31(2,:); 

kg31(2,:)=templ; 

temp2=kg21(:,l); 

kg31(:,l)=kg  31(:,2); 

kg31(:,2)=temp2; 

pg3i=pgi; 

temp=pg31(l); 

Pg3 1(1  )=pg3 1  (2) ; 

Pg31(2)=temp; 

m3c=mg31(dpn:dpn*(disp-l)*num,dpn:dpn*(disp-l)*num); 
k3c=kg31(dpn:dpn*(disp-l)*num,dpn:dpn*(disp-l)*num); 
p3c=pg31(dpn:dpn*(disp-l)*num); 
c3c  =  alpha*k2c  +  beta*m2c; 
end 
end 
end 

%Static  solution 

w01(dpn:dpn*(disp-l)*num)=k0c\p0c; 
if  order  >  0 

wl  I(dpn:dpn*(disp-l)*num)=-k0c\(klc*w01(dpn:dpn*(disp-l)*num)-plc); 
if  order  >  1 

w21(dpn:dpn*(disp-l)*num)=-k0c\(klc*wll(dpn:dpn*(disp-l)*num)+k2c*w01(dpn:dpn*(disp- 
l)*num)-p2c); 
if  order  >  2 

w31(dpn:dpn*(disp-l)*num)=-k0c\(klc*w21(dpn:dpn*(disp-l)*num)+k2c*wll(dpn:dpn*(disp- 

I)*num)+k3c*w01(dpn:dpn*(disp-l)*num)-p3c); 

end 

end 

end 

%  Resort  displacements 

w0=w01; 

temp=wO(l); 

w0(l)=w0(2); 

w0(2)=temp; 

if  order  >  0 

wl=wll; 

temp=wl(l); 

wl(l)=wl(2); 
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wl(2)=temp; 

if  order  >  1 

w2=w21; 
temp=w2(l); 
w2(l)=w2(2); 
w2(2)=temp; 
if  order  >  2 

w3=w31; 

temp=w3(l); 

w3(l)=w3(2); 

w3(2)=temp; 

end 

end 

end 

start=els(l).rl; 
for  i=l:num 

nodeids=[(disp-l)*i-(disp-2),  (disp-l)*i+(disp-2),  (disp-l)*i]; 

dO=GetDisplacements(wO,dpn,nodeids); 

if  order  >  0 

dl=GetDisplacements(wl,dpn,nodeids); 
if  order  >  f 

d2=GetDisplacements(w2,dpn,nodeids); 
if  order  >  2 

d3=GetDisplacements(w3,dpn,nodeids); 

end 

end 

end 

for  j=0:chunk 

%  ds=(diag([els(i).len  1  els(i).len  1  els(i).len  l])/els(i).eta)*dO; 
x((i-l)*chunk+j+l)=start  +  (j/chunk)*els(i).len; 

ws((i-f)*chunk+j+l,f)=memshape(dO,x((i-l)*chunk+j+l),els(i).eps,els(i).r,els(i).len); 
if  order  >  0 

ws((i-l)*chunk+j+l,2)=memshape(d0+dl,x((i-l)*chunk+j+l),els(i).eps,els(i).r,els(i).len); 
if  order  >  1 

ws((i-i)*chunk+j+l,3)=memshape(d0+dl+d2,x((i-l)*chunk+j+l),els(i).eps,els(i).r,els(i).len); 

if  order  >  2 

ws((i-l)*chunk+j+l,4)=memshape(d0+dl+d2+d3,x((i- 

l)*chunk+j+l),els(i).eps,els(i).r,els(i).len); 

end 

end 

end 

end 

start=start+els(i).len; 


if  dyn  ==  0 
if  verbose  >  0 
figure 

wslamda=ws/(633*10A-9); 
plot(x,wslamda); 
title('Static  Shape'); 
legend('wO',  'wl',  'w2'); 
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figure 

wslamda2(:,order+l)=(ws(:,order+l))/(633*10A-9); 
plot(x,wslamda2(:  ,order+ 1 )); 
title('Static  Shape'); 

%  legend('w'); 

%  legend('FEM  Solution'); 
end 

xr  =  x; 
vr  =  ws; 
dr  =  0; 
return 
else 

%Dynamic 

v01=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 
d01=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 
if  order  >  0 

vl  l=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 
dl  l=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 
if  order  >  1 

v2=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 
d2=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 
if  order  >  2 

v3=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 

d3=zeros(dpn*((disp-l)*num+l),dpn*((disp-l)*num+l)); 

end 

end 

end 


if  mode  ==  0 

[v01(dpn:dpn*(disp-l)*num,dpn:dpn*(disp-l)*num),d01(dpn:dpn*(disp-l)*num,dpn:dpn*(disp- 

I)*num)]=eig(k0c,m0c); 

%  Resort  displacements 
v0=v01; 
templ=v0(l,:); 
v0(l,:)=v0(2,:); 
v0(2,:)=templ; 

[vc0,dc0]=sorteigs(v0,d0,0,mx); 

if  order>0 
for  i=l:length(dcO) 

omegal(i)=(vO(:,i)'*(kgl-dcO(i)*mgl)*vO(:,i)); 
den=v0(:,i)'*mg0*v0(:,i); 
if  den  ==  0 
omegal(i)=0; 
else 

omega  1  (i)=omega  1  (i)/2/dc0(i)/den; 
end 
end 

if  order>l 
for  i=l:length(dcO) 

omega2(i)=(vO(:,i)'*(kg2-dcO(i)*(omegal(i)A2*mgO+2*omegal(i)*mgl+mg2))*vO(:,i)); 

den=v0(:,i)'*mg0*v0(:,i); 
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if  den  ==  0 
omega2(i)=0; 
else 

omega2(i)=omega2(i)/2/dc0(i)/den; 

end 

end 

if  order>2 
for  i=l:length(dcO) 
omega3(i)=(v0(:,i)'*(kg3- 

dc0(i)*(2*omegal(i)*omega2(i)*mg0+(2*omega2(i)+omegal(i)A2)*mgl+2*omegal(i)*mg2+mg3))*v0(:,i 

)); 

den=vO(:,i)'*mgO*vO(:,i); 

if  den  ==  0 
omega3(i)=0; 
else 

omega3(i)=omega3(i)/2/dc0(i)/den; 

end 

end 

end 

end 

end 


if  verbose>0 

if  length(dc0)>0 
figure; 

[x,w]=plotmemshape(vc0,dpn,els,2); 
title('Clamped  VO  Modes'); 
figure; 
bar(dcO); 

title( 'Clamped  DO  -  Eigenvalues'); 
figure; 

bar(sqrt(dcO/dcO(l))); 

title) 'Clamped  First  Order  Frequencies  (Scaled)'); 
end 

if  order  >  0 
figure; 

bar(omegal); 

title)' 1st  Order  Frequency  Corrections'); 

if  order  >  1 
figure; 

bar(omega2); 

title('2nd  Order  Frequency  Corrections'); 

if  order  >  2 
figure; 

bar(omega3); 

title('2nd  Order  Frequency  Corrections'); 
end 
end 
end 
else 

start =e  Is)  l).rl; 
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for  j=l:num 
for  k=0:chunk 

x((j-l)*chunk+k+l)=start+(k/chunk)*els(j).len; 
for  i=l:length(dcO) 

ve=GetDisplacements(vcO(:,i),  dpn,  [(j-l)*2+l  j*2+l  j*2]); 

w((j-l)*chunk+k+l  ,i)=memshape(ve,x((j-l)*chunk+k+l),els(j).eps,els(j).r,els(j).len); 
end 

end 

start=start+els(j).len; 

end 

end 
vr  =  w; 
xr  =  x; 
dr  =  dcO; 
else 

%Force  Response 
[vO  1  ,d0 1  ]=eig(k0c,m0c); 

[vOc,dOc]=sorteigs(vO  1  ,d0  l,0,mx); 

for  k=l :  1 : 100 1 
delta(k)=0.00001*(k-501); 

a20=v0c(:,l)'*plc/(v0c(:,l)'*(klc+sqrt(-d0c(l))*clc-d0c(l)*mlc- 
2*delta(k)*sqrt(d0c(l))*m0c)*v0c(:,l)); 
if  order  >  0 

bl(:,k)=plc-a20*(klc+sqrt(-d0c(l))*clc-d0c(l)*mlc-2*delta(k)*sqrt(d0c(l))*m0c)*v0c(:,l); 

av=(k0c-d0c(mode)*m0c)*v0c; 

c(:,k)=pinv(av)*bl(:,k); 

c(l,k)=a20; 

end 

end 

if  verbose  >  0 
figure; 

plot(delta,c); 

xlabel('\delta'); 

end 
vr  =  c; 
xr  =  delta; 
dr  =  dOc; 
end 
end 
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Another  file  that  performs  some  calculations  and  creates  plots  for  the  Zernike  coefficients 
is  included  here.  The  file  was  reduced  for  use  in  the  final  code  and  renamed 
runZernsingle. 

function  runZern(num,verb,config) 
if  verb  >  0 

%  Calculations  for  no  pressure  and  changing  voltage  and  constant  tension 
volts= [0:50:1 000] ; 
for  i=l:length( volts) 

[v,x,temp,eps]=testr2(num,0,1000,volts(i),2,0,0,0,0,0,config); 
a{  1  }.coef(i,:)=GetAxiZerntest(x',2*v(:,3)'); 
a{  1  }.eps(i,:)=eps 
end 

%  Define  voltage  vector  for  output 
Voltage  =  volts'; 

%  Calculations  for  .01  pressure  and  changing  voltage  and  constant  tension 
for  i=  1  :length(  volts) 

[v,x, temp, eps]=testr2(num, 0.01, 1000, volts(i), 2, 0,0,0,0,0,config); 
a{  2 }  .coef(i,:)=GetAxiZerntest(x',2*v(:,3)'); 
a{2}.eps(i,:)=eps; 
end 

%  Graphing  changing  voltage  for  above  cases 
if  verb  >  0 
figure 

subplot(2,l,l) 

plot(volts,a{  1 }. coef(:,l),'k--', volts, a{  l}.coef(:,2),’k-.',volts,a{  1 } .coef(:,3),'k-r); 

a{  1  }.title=’Pressure  =  O'; 

title('Pressure  =  O'); 

legend([’Z_l  ';  'Z_5  ';  'Z_{  13}’]); 

subplot(2,l,2) 

plot(volts,a{2}.coef(:,l),'k--',volts,a{2}.coef(:,2),'k-.',volts,a{2}.coef(:,3),'k-'); 
a{2}.title=’Pressure  =  0.01'; 
title('Pressure  =  0.01'); 
legend([’Z_l  ';  ’Z_5  ';  'Z_{  13}']); 
xlabel('Voltage  Deviations  (T=1000)'); 
end 

%  Calculations  for  no  voltage,  low  pressure  and  changing  tension 
tens=[  1000: 1000:7000]; 
for  i=l:length(tens) 

[v,x,temp,eps]=testr2(num,0. 01, tens(i), 0,2, 0,0,0, 0,0,config); 
a{3}.coef(i,:)=GetAxiZemtest(x',2*v(:,3)'); 
a{3}.eps(i,:)=eps; 
end 

%  Calculations  for  .1  pressure  and  100  voltage  and  changing  tension 
for  i=l:length(tens) 

[v,x,temp,eps]=testr2(num,0.1,tens(i),  100,2, 0,0, 0,0,0,config); 
a{4}.coef(i,:)=GetAxiZerntest(x',2*v(:,3)'); 
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a{4}.eps(i,:)=eps; 

end 

%  Graphing  changing  tensions  for  above  cases 
if  verb  >  0 
figure 

subplot(2,f  ,f) 

plot(tens,a{3}.coef(:,l),'k— ',tens,a{3}.coef(:,2),'k-.',tens,a{3}.coef(:,3),'k-'); 

a{  3 }  .title='d_{  3 1 }  V  =  0'; 

title('d_{  3 1 }  V  =  O'); 

legend(['Z_f  'Z_5  'Z_{  f3}’]); 

subplot(2,f  ,2) 

plot(tens,a{4}.coef(:,l),'k— ',tens,a{4}.coef(:,2),'k-.',tens,a{4}.coef(:,3),'k-'); 
a{4}.title=’d_{31}V  =  fOA{-6}’; 
title('d_{31}V  =  10A{-6}’); 
legend(['Z_l  'Z_5  ’Z_{  13}’]); 

xlabelfTension  Deviations  (P=0.0 1)’); 
end 


%  Calculations  for  changing  pressure  and  no  voltage  and  1000  tension 
pres=[0:0. 001:0.01]; 
for  i=l:length(pres) 

[v,x, temp, eps]=testr2(num,pres(i),  1000,0, 2, 0,0,0,0,0,config); 
a{5  }.coef(i,:)=GetAxiZerntest(x',2*v(:,3)'); 
a{5  }.eps(i,:)=eps; 
end 

%  Calculations  for  changing  pressure  and  100  voltage  and  1000  tension 
for  i=l:length(pres) 

[v,x,temp,eps]=testr2(num,pres(i),  1000, 100,2, 0,0,0,0,0,config); 
a{  6 }  .coef(i,:)=GetAxiZerntest(x',2*v(:,3)'); 
a{6}.eps(i,:)=eps; 
end 

%  Graphing  changing  pressure  for  above  cases 
if  verb  >  0 
figure 

subplot(2,l,l) 

plot(pres,a{5  }.coef(:,l),’k— ',pres,a{5}.coef(:,2),'k-.',pres,a{5  }.coef(:,3),'k-'); 

a{  5 }  .title='d_{  3 1 }  V  =  O'; 

title(’d_{  3 1 }  V  =  O'); 

legend] [’Z_l  ';  'Z_5  ’;  'Z_{  13}’]); 

subplot(2,l,2) 

plot(pres,a{6}.coef(:,l),'k— ',pres,a{6}.coef(:,2),'k-.',pres,a{6}.coef(:,3),'k-'); 
a{6}. title  =  'd_{  3 1 }  V  =  10A  { -6}'; 
title]  'd_  { 3 1 }  V  =  10A{-6}'); 
legend] [’Z_l  ';  'Z_5  ’;  ’Z_{  13}’]); 
xlabel('Pressure  Deviations  (T=1000)'); 
end 
end 
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The  rest  of  the  files  are  ones  created  to  make  the  program  run  smoother  and  complete  all 
the  necessary  actions  by  running  only  one  program.  That  file  is  included  first,  followed 
by  other  files  that  are  used  for  calculations  to  create  the  7  Figures. 

function  Program 
%  Clearing 
clc 
clear; 

%  Assignment  for  while  loops  and  other  functions 

global  zfe 

verbl  =  ’q'; 

verb2  =  ’q'; 

zerocond  =  'q'; 

rover  =  [0;.03:3]'; 

rmin  =  rover/max(rover); 

ang  =  [-pi/2:pi/12:3*pi/2]'; 

mkpts  =  1;  %  Zero  to  Stop  Experimental  Plots,  1  to  see  them 

verb  =  0;  %  Change  to  1  (one)  if  you  want  to  see  the  plots  created  by  original  FEM  code 

%  Questions  that  need  to  be  answered  to  get  program  to  run 
%  nurnl  =  input(’How  many  symmetric  elements  do  you  want?  '); 

%  config  =  input('What  configuration  do  you  want  to  activate(0  or  1  for  even  elements,  2-5  for  odd 
elements --5  is  1"  circle  in  center)?  '); 

%  volt  =  input('What  voltage  do  you  want  to  apply  to  the  membrane(for  comparison  -600  to  600  in  300V 
increments)? '); 

%  %  p  =  input('What  pressure  do  you  want  to  apply  to  the  membrane?  '); 

%  t  =  input('What  tension  is  placed  on  the  membrane  (use  1000)?  '); 

%  Setting  values  so  questions  don't  have  to  be  there— Comment  if  you  have  questions  active 
numl  =  9; 
config  =5 ; 
volt=-600; 

p=0; 

t=80000000; 

%  Stuff  for  graphs 
Element  =  num2str(numl); 

V=num2str(volt); 

P=num2str(p); 

N=num2str(t); 

%  Different  solutions  for  the  zero  voltage  case 
if  (volt  ==  0) 

while  (zerocond  ~=  'b')  &  (zerocond  ~=  'm')  &  (zerocond  ~=  'a') 

zerocond  =  input('Which  0  V  condition  do  you  want  (before,  mid,  after)?  ','s'); 
if  (zerocond  ==  'b') 
zerocon  =  1 ; 
elseif  (zerocond  ==  ’m') 
zerocon  =  2; 
elseif  (zerocond  ==  'a') 
zerocon  =  3 ; 
else 
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'Please  answer  with  proper  input' 
end 
end 
else 

zerocon  =  0; 
end 


%  FEM  output 
runZern(num  1 ,  verb.config) 

Zcoeff  =  runzernsingle(numl,verb,config,volt,p,t); 
for  i=l:6 

zfe(i)  =  Zcoeff(l,i); 
end 

testr2(numl,p,t,volt,2,verb,0,0,0,0,config); 

[v,x,d,eps,youngsperlayer,thickperlayer]  =  testr2(numl,p,t,volt,2,verb,0,0,0,0,config); 
Emem  =  num2str(youngsperlayer(l,l)); 
tmem  =  num2str(thickperlayer(l,l)); 

%  Experimental  output  and  deflection  vectors 
if  volt  ==  0 
if  zerocon  ==  1 
zerocenter 
elseif  zerocon  ==  2 
zerocentermid 
elseif  zerocon  ==  3 
zerocenterafter 
end 

elseif  volt  ==  300 
threeOOcenter 
elseif  volt  ==  -300 
neg300center 
elseif  volt  ==  -600 
neg600center 
elseif  volt  ==  600 
sixhundredcenter 

end 

%  Running  outside  script  for  deflection  calculations 

[wfem,wexp,wexpall,wexpnon]  =  deflection(zfe,z,rover,rmin,ang,volt,mkpts,t); 

wexpallmod  =  wexpall; 
wexpnonmod  =  wexpnon; 

%  Deflection  plots 
wvs  =  wexp-wfem; 
angl  =  -pi/2; 
for  i  =  l:size(rmin) 
wexpallred(i)  =  0; 
wexpnonred(i)  =  0; 
for  j  =  l:size(ang) 
if  (angl  >=  pi/2) 
for  k  =  l:size(rmin) 
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for  m  =  l:size(ang) 
test  =  abs(ang(m)  -  angl); 
if  (test  <  pi/24) 
n  =  m; 
end 
end 

wexpallmod(k,n)  =  -wexpallmod(k,n); 
wexpnonmod(k,n)  =  -wexpnonmod(k.n); 
end 
end 

angl  =  angl  +  pi/12; 

wexpallred(i)  =  wexpallred(i)  +  wexpallmod(i,j); 
wexpnonred(i)  =  wexpnonred(i)  +  wexpnonmod(i.j); 
end 
end 

wexpallred'; 

wexpnonred'; 

figure 

subplot(2,l,l),  plot(rover,wexp,rover,wfem),title('Experimental  and 

FEM’, 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'b'),xlabel('Radial  dimension  edge  to  center'),... 
ylabel(’Wavefront  deflection') 

subplot(2,l,2),  plot(rover,wvs),title(’Difference  between  Experimental  and 

FEM’, FontSize', 12, 'FontWeight', 'bold', 'Color', 'g'),xlabel('Radial  dimension  edge  to  center'),... 

ylabel(’Wavefront  deflection') 
figure 

subplot(2,l,l),plot(rover,wfem),title(['Finite  Element  Solution  Voltage  =  ' 

V],'FontSize', 12, 'FontWeight', 'bold', 'Color', 'r'),xlabel('Radial  dimension  edge  to  center'),... 
ylabel(’Wavefront  deflection') 

subplot(2,l,2),plot(rover,wexp),title('Experimental  Reduced  (Symmetric) 

Solution', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'r'),xlabel('Radial  dimension  edge  to  center'),... 

ylabel('Wavefront  deflection') 
figure 

subplot(2,l,l),plot(rover,wexpall),title('All  35  Zernike  from  -Pi/2  to 

3*Pi/2', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'y'),xlabel('Radial  dimension  edge  to  center'),... 
ylabel('Wavefront  deflection') 

subplot(2,l,2),plot(rover,wexpnon),title('30  Non-Axisymmetric  Zernike 

Coefficents', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'b'),xlabel('Radial  dimension  edge  to  center'),... 

ylabel('Wavefront  deflection') 
figure 

subplot(2,l,l),plot(rover,wexpallred),title(’Reduced 

35', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'c'),,xlabel('Radial  dimension  edge  to  center'),... 

ylabel('Wavefront  deflection') 
subplot(2,l,2),plot(rover,  wexpnonred), title(’Reduced 

30', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'm'),xlabel('Radial  dimension  edge  to  center'),... 

ylabel('Wavefront  deflection') 
figure 

subplot(2,l,l),plot(rover,wexpallmod),title(’Modified  35 

Zernike', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'r'),xlabel('Radial  dimension  edge  to  center'),... 

ylabel(  Wavefront  deflection') 
subplot(2,l,2),plot(rover,wexpnonmod),title(’Modified  30 

Zernike', 'FontSize',  12, 'FontWeight', 'bold', 'Color', 'b'),xlabel(’Radial  dimension  edge  to  center'),... 
ylabel('Wavefront  deflection') 
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This  file  is  a  sub  file  to  the  above  Program  file  and  is  used  to  calculate  the  deflections 
based  on  the  Zemike  polynomials. 

function  [wfem,wexp,wexpall,wexpnon]  =  deflection(zfe,z,rover,rmin,ang,volt,mkpts,tension) 

rover  =  [0:.03:3]'; 

rmin  =  rover/max(rover); 

ang  =  [-pi/2:pi/12:3*pi/2]'; 

%  Code  for  deflections  along  any  center  line  (FEM  and  experiment--  5  zernike  numbers) 
for  i  =  l:size(rmin) 

wfem(i,:)  =  zfe(2)*(2*((rover(i,:))A2)-l)  +  zfe(3)*(6*((rover(i,:))A4)-6*((rover(i,:))A2)+l)  + 
zfe(4)*(20*((rover(i,:))A6)-30*((rover(i,:))A4)+12*((rover(i,:))A2)-l)  +  ... 

zfe(5)*(70*((rover(i,:))A8)-140*((rover(i,:))A6)+90*((rover(i,:))A4)-20*((rover(i,:))A2)+l)  +  ... 
zfe(6)*(252*(rover(i,:))A10-630*(rover(i,:))A8+560*(rover(i,:))A6- 
210*(rover(i,:))A4+30*(rover(i,:))A2-l); 

wexp(i,:)  =  z(3,:)*(2*(rmin(i,:))A2-l)  +  z(8,:)*(6*(rmin(i,:))A4-6*(rmin(i,:))A2+l)  + 
z(15,:)*(20*(rmin(i,:))A6-30*(rmin(i,:))A4+12*(rmin(i,:))A2-l)  +  ... 

z(24,:)*(70*(rmin(i,:))A8-140*(rmin(i,:))A6+90*(rmin(i,:))A4-20*(rmin(i,:))A2+l)  +  ... 
z(35,:)*(252*(rmin(i,:))A10-630*(rmin(i,:))A8+560*(rmin(i,:))A6-210*(rmin(i,:))A4+30*(rmin(i,:))A2- 

1); 

end 

if  (volt  <  0) 

wexp  =  wexp  *  -1; 
end 

%  Code  for  Experimental  deflection  for  angles  from  0  to  Pi  (35  zernike  numbers) 
if  volt  ==  0 
if  zerocon  ==  1 

[z]=M2_Tlb_0v(mkpts); 
elseif  zerocon  ==  2 
M2_T  1  e_0V  (mkpts) ; 
elseif  zerocon  ==  3 
M2_T  lh_0  V  ( mkpts); 
end 

elseif  volt  ==  300 

[z]  =M2_T  1  c_300v(mkpts) ; 
elseif  volt  ==  -300 

[z]  =M2_T  1  f_Ne  g300v(mkpts) ; 
elseif  volt  ==  -600 

[z]=M2_Tlg_Neg600v(mkpts); 
elseif  volt  ==  600 

[z]  =M2_T  1  d_600v(mkpts) ; 

end 

for  i  =  1:  size  (rmin) 
for  j  =  l:size(ang) 

wexpall(i.j)  =  z(l,:)*(rmin(i)*cos(ang(j)))  +  z(2,:)*(rmin(i)*sin(ang(j)))  +  z(3,:)*(2*rmin(i)A2-l)  + 
z(4,:)*(rmin(i)A2*2*cos(2*ang(j)))  +  ... 

z(5,:)*(rmin(i)A2*2*sin(2*ang(j)))  +  z(6,:)*((3*rmin(i)A2-2)*rmin(i)*cos(ang(j)))  + 
z(7,:)*((3*rmin(i)A2-2)*rmin(i)*sin(ang(j)))  +  ... 

z(8,:)*(6*(rmin(i))A4-6*(rmin(i))A2+l)  +  z(9,:)*(rmin(i)A3*cos(3*ang(j)))  + 
z(10,:)*(rmin(i)A3*sin(3*ang(j)))  +  ... 

z(l  l,:)*((4*rmin(i)A2-3)*rmin(i)A2*cos(2*ang(j)))  +  z(12,:)*((4*rmin(i)A2- 
3)*rmin(i)A2*sin(2*ang(j)))  +  ... 
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z(i3,:)*((iO*rmin(i)A4-f2*rmin(i)A2+3)*rmin(i)*cos(ang(j)))  +  z(14,:)*((10*rmin(i)A4- 
f2*rmin(i)A2+3)*rmin(i)*sin(ang(j)))  +  ... 

z(15,:)*(20*(rmin(i))A6-30*(rmin(i))A4+12*(rmin(i))A2-l)  +  z(16,:)*(rmin(i)A4*cos(4*ang(j)))  + 
z(17,:)*(rmin(i)A4*sin(4*ang(j)))  +  ... 

z(18,:)*((5*rmin(i)A2-4)*rmin(i)A3*cos(3*ang(j)))  +  z(19,:)*((5*rmin(i)A2- 
4)*rmin(i)A3*sin(3*ang(j)))  +  ... 

z(20,:)*((f5*rmin(i)A4-20*rmin(i)A2+6)*rmin(i)A2*cos(2*ang(j)))  +  z(21,:)*((15*rmin(i)A4- 
20*rmin(i)A2+6)*rmin(i)A2*sin(2*ang(j)))  +  ... 

z(22,:)*((35*rmin(i)A6-60*rmin(i)A4+30*rmin(i)A2-4)*rmin(i)*cos(ang(j)))  + 
z(23,:)*((35*rmin(i)A6-60*rmin(i)A4+30*rmin(i)A2-4)*rmin(i)*sin(ang(j)))  +  ... 

z(24,:)*((70*(rmin(i))A8-140*(rmin(i))A6+90*(rmin(i))A4-20*(rmin(i))A2+l))  + 
z(25,:)*(rmin(i)A5*cos(5*ang(j)))  +  ... 

z(26,:)*(rmin(i)A5*sin(5*ang(j)))  +  z(27,:)*((6*rmin(i)A2-5)*rmin(i)A4*cos(4*ang(j)))  +  ... 
z(28,:)*((6*rmin(i)A2-5)*rmin(i)A4*sin(4*ang(j)))  +  z(29,:)*((21*rmin(i)A4- 
30*rmin(i)A2+10)*rmin(i)A3*cos(3*ang(j)))  +  ... 

z(30,:)*((21*rmin(i)A4-30*rmin(i)A2+10)*rmin(i)A3*sin(3*ang(j)))  +  z(31,:)*((56*rmin(i)A6- 
105*rmin(i)A4+60*rmin(i)A2-10)*rmin(i)A2*cos(2*ang(j)))  +  ... 

z(32,:)*((56*rmin(i)A6-105*rmin(i)A4+60*rmin(i)A2-10)*rmin(i)A2*sin(2*ang(j)))  +  ... 
z(33,:)*((126*rmin(i)A8-280*rmin(i)A6+210*rmin(i)A4-60*rmin(i)A2+5)*rmin(i)*cos(ang(j)))  +  ... 
z(34,:)*((126*rmin(i)A8-280*rmin(i)A6+210*rmin(i)A4-60*rmin(i)A2+5)*rmin(i)*sin(ang(j)))  +  ... 
z(35,:)*(252*rmin(i)A10-630*rmin(i)A8+560*rmin(i)A6-210*rmin(i)A4+30*rmin(i)A2-l); 

end 

end 

%  3 1  Non  symmetric  zernike  numbers 
for  i  =  l:size(rmin) 
for  j  =  f:size(ang) 

wexpnon(i.j)  =  z(f  ,:)*(rmin(i)*cos(ang(j)))  +  z(2,:)*(rmin(i)*sin(ang(j)))  + 
z(4,:)*(rmin(i)A2*2*cos(2*ang(j)))  +  ... 

z(5,:)*(rmin(i)A2*2*sin(2*ang(j)))  +  z(6,:)*((3*rmin(i)A2-2)*rmin(i)*cos(ang(j)))  + 
z(7,:)*((3*rmin(i)A2-2)*rmin(i)*sin(ang(j)))  +  ... 

z(9,:)*(rmin(i)A3*cos(3*ang(j)))  +  z(fO,:)*(rmin(i)A3*sin(3*ang(j)))  +  ... 
z(f  f  ,:)*((4*rmin(i)A2-3)*rmin(i)A2*cos(2*ang(j)))  +  z(f2,:)*((4*rmin(i)A2- 

3) *rmin(i)A2*sin(2*ang(j)))  +  ... 

z(i3,:)*((iO*rmm(i)A4-f2*rmin(i)A2+3)*rmin(i)*cos(ang(j)))  +  z(f4,:)*((fO*rmin(i)A4- 
f2*rmin(i)A2+3)*rmin  (i)*sin(ang(j)))  +  ... 

z(f6,:)*(rmin(i)A4*cos(4*ang(j)))  +  z(i7,:)*(rmin(i)A4*sin(4*ang(j)))  +  ... 
z(f8,:)*((5*rmin(i)A2-4)*rmin(i)A3*cos(3*ang(j)))  +  z(i9,:)*((5*rmin(i)A2- 

4) *rmin(i)A3*sin(3*ang(j)))  +  ... 

z(20,:)*((f5*rmin(i)A4-20*rmin(i)A2+6)*rmin(i)A2*cos(2*ang(j)))  +  z(2f,:)*((i5*rmin(i)A4- 
20*rmin(i)A2+6)*rmin(i)A2*sin(2*ang(j)))  +  ... 

z(22,:)*((35*rmm(i)A6-60*rmin(i)A4+30*rmin(i)A2-4)*rmin(i)*cos(ang(j)))  + 
z(23,:)*((35*rmin(i)A6-60*rmin(i)A4+30*rmin(i)A2-4)*rmin(i)*sin(ang(j)))  +  ... 

z(25,:)*(rmin(i)A5*cos(5*ang(j)))  +  z(26,:)*(rmin(i)A5*sin(5*ang(j)))  +  z(27,:)*((6*rmin(i)A2- 

5) *rmin(i)A4*cos(4*ang(j)))  + ... 

z(28,:)*((6*rmin(i)A2-5)*rmin(i)A4*sin(4*ang(j)))  +  z(29,:)*((2f  *rmin(i)A4- 
30*rmin(i)A2+i0)*rmin(i)A3*cos(3*ang(j)))  +  ... 

z(30,:)*((2f*rmin(i)A4-30*rmin(i)A2+f0)*rmin(i)A3*sin(3*ang(j)))  +  z(3i,:)*((56*rmin(i)A6- 
f05*rmin(i)A4+60*rmin(i)A2-f0)*rmin(i)A2*cos(2*ang(j)))  +  ... 

z(32,:)*((56*rmin(i)A6-f05*rmin(i)A4+60*rmin(i)A2-f0)*rmin(i)A2*sin(2*ang(j)))  +  ... 
z(33,:)*((f26*rmin(i)A8-280*rmin(i)A6+2i0*rmin(i)A4-60*rmin(i)A2+5)*rmin(i)*cos(ang(j)))  +  ... 
z(34,:)*((f26*rmin(i)A8-280*rmin(i)A6+2f0*rmin(i)A4-60*rmin(i)A2+5)*rmin(i)*sin(ang(j))); 

end 

end 
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Appendix  B:  Code  Operation 


For  most  users  of  the  finite  element  code,  the  function  file  above  called  Program 
was  created  that  will  ask  simple  questions  about  the  system.  The  user  answers  these 
questions  and  the  program  runs  for  about  30  seconds  while  calculations  are  completed 
and  various  plots  are  produced.  The  program  will  currently  produce  the  following  plots: 
Figure  1 — 4  plots  of  Original  Experimental  Data  (Zemike  and  Raw  data) 

Figure  2 — 4  plots  of  Reduced  Experimental  Data 

Figure  3 — Experimental  and  FEM  solution  on  same  plot.  Difference  of  two  solutions 
Figure  4 — Experimental  and  FEM  solution  on  different  plot 

Figure  5 — All  35  Experimental  Zemike  numbers  along  different  center  lines  (different 
angles) 

Figure  6 — Same  as  Figure  5  for  just  the  non-axisymmetric  Zemike  numbers 
Figure  7 — The  addition  of  the  many  centerlines  for  all  angles  around  the  membrane 
This  concludes  the  basic  guide  for  the  everyday  user  of  the  program. 

For  a  more  advanced  user,  it  may  be  necessary  to  understand  the  code  in  more 
depth.  The  files  included  in  Appendix  A  are  the  main  files  that  a  user  may  need  to  make 
changes  to  for  specific  conditions.  All  the  necessary  files  (as  well  as  other  files  created 
during  the  development  of  the  original  code)  are  included  on  a  CD  with  this  thesis. 

Some  of  these  files  will  be  discussed  now  to  supply  the  reader  with  more 
understanding  and  knowledge.  Basic  understanding  of  the  input  commands  in  MATFAB 
is  assumed. 

As  stated  before,  most  of  the  important  code  is  included  in  Appendix  A.  The  first  file 
that  will  be  discussed  is  testr2.  This  file  contains  the  material  properties  and  has  inputs 
related  to  the  configuration  and  number  of  elements.  The  size  of  the  radius  is  defined  at 
the  top.  It  is  to  be  noted  that  the  first  number  cannot  be  set  to  zero  or  there  will  be  an 
overflow  in  the  calculations.  The  layup  defined  is  the  polymer  membrane  followed  by 
the  piezo  layer;  therefore,  the  material  properties  are  listed  in  this  fashion.  For  example, 
there  is  one  line  that  reads: 

youngsperlayer=diag(  [  1 00  261000]); 

The  first  number  in  the  brackets  is  the  modulus  for  the  membrane  and  the  second 
is  for  the  piezo  layer.  This  is  the  same  for  all  the  material  properties.  Mostly,  only  the 
modulus,  thickness  and  possibly  the  density  (if  running  a  dynamic  solution)  should  be 
changed  in  this  file.  The  input  variables  are  listed  in  the  top  function  line: 

testr2(numel,p,t,v,order,verb,dyn,mode,alpha,beta,config); 

They  are  in  order: 

Numel:  number  of  elements  along  one  radial  line 
P:  Pressure  term  (usually  this  should  be  assigned  zero) 
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T:  Tension  term 

V:  Voltage  applied  to  the  region 

Order:  Perturbation  order  desired  (this  should  usually  be  set  to  2) 

Verb:  Verbose  command  If  set  to  anything  but  zero,  the  program  has  plots  that  it  will 
output  that  were  used  in  the  dissertation. 

Dyn:  Set  to  1  if  a  dynamic  solution  is  desired 

Mode:  Used  if  a  dynamic  solution  is  being  solved.  Wasn’t  used  during  this  study. 

Alpha,  Beta:  Not  used  in  static  study. 

Config:  Allows  the  user  to  quickly  pick  a  different  configuration  of  the  piezo  activation 
zones.  For  example,  configuration  5  is  approximately  a  1”  center  circle  (for  the  plots 
since  will  only  represent  half  of  this  on  the  plots  due  to  the  modeling  of  a  radial  line).  As 
it  is  noticed  near  the  bottom  of  this  file,  the  configuration  is  defined  by: 

elseif  config  ==  5 

volts=v*[ones(2,numel/3)  zeros(2,numel/3)  zeros(2,numel/3)  zeros(2,numel/3)  zeros(2,numel/3)]; 

Notice  there  are  five  zones  over  the  membrane  with  a  3”  radius.  Therefore  each  zone  is 
3/5”  in  length.  Most  of  the  zones  are  zero  except  for  the  most  left  one  which  represents 
the  center.  If  another  one  of  these  was  changed  to  ones  (let’s  say  the  third  one),  there 
would  be  a  6/5”  circle  in  the  center  and  a  ring  with  an  inner  diameter  of  6/5”  and  outer 
diameter  of  9/5”  that  would  be  activated.  With  enough  imagination  and  slight  changes  in 
the  code  can  create  all  sorts  of  symmetric  activation  zones.  It  is  advisable  to  keep  the 
zones  at  a  reasonable  size,  not  allowing  them  to  get  too  small,  just  to  prevent  losing 
information  across  each  zone.  Usually,  an  odd  number  of  elements  is  best  and  depending 
on  whether  odd  or  even  number  are  used  will  confine  to  which  configuration  can  be  used. 

As  with  all  the  original  finite  element  codes,  it  is  best  to  run  the  function  script 
Program  rather  than  running  one  at  a  time.  The  other  program  that  is  needed  is 
runzernsingle.  It  will  actually  run  the  above  file  ( testr2 ),  so  it  requires  inputs  of  the 
number  of  elements,  verbose  command  (set  to  zero)  and  the  desired  configuration.  There 
is  nothing  to  be  changed  in  this  file. 

The  Program  and  Deflection  files  are  the  main  front  end  files.  The  deflection 
script  will  probably  never  need  too  much  modification.  The  program  script  runs  the 
show  and  it  can  be  setup  to  just  run  or  have  questions  setup  on  the  command  line 
allowing  the  user  to  answer  the  questions.  All  the  definitions  of  the  properties  are  near 
the  top  of  the  file.  This  program  can  be  run  from  the  command  line  by  typing  program  if 
the  directory  is  set  properly.  Otherwise,  there  is  a  run  button  if  the  file  is  open  in  Matlab. 
There  is  not  much  to  this  file  and  modifications  to  the  graphs  can  be  completed  at  the 
bottom.  The  center  area  is  running  the  other  scripts  defined  above.  Any  basic  operation 
of  the  program  should  be  directed  to  Capt  Brian  Lutz. 
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